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FOREWORD 

The objective of Ocean Floor Rock Embedment Mechanics (OFREM) research is to 
develop a better understanding of interaction of rock media and ballistically embedded projectiles 
that ultimately can serve as ocean floor anchors. Current tasks are to: (a) determine the model 
scale effects for penetration into intact and jointed rock, (b) mathematically simulate the 
penetration process, (c) determine projectile extraction resistance and projectile/target interaction, 
and (d) mathematically simulate the projectile extraction process. 

From FY88 through FY90, the University of California at Berkeley under contract to 
NCEL was tasked to formulate analytical models for predicting projectile anchor extraction in 
rock masses. Due to personnel turnover, the originally proposed analytical work based upon the 
use of the DYNA codes from Lawrence Livermore National Laboratory was discontinued. At 
about the same time, the DYNA codes were imported for use on unrelated projects in the 6.2 
NAVFAC/NCEL Exploratory Development program. The present project resumes the originally 
proposed research based upon these codes, addressing the problems of numerical simulation of 
the penetration process, projectile/target interaction, and prediction of extraction resistance. 

A similar project report for this study was submitted by the first author in partial 
fulfillment of the requirements for the degree of Master of Science in Mechanical Engineering 
at California State University, Northridge in May 1992; the title of the thesis is "Stress 
Surrounding Projectile Anchors." 
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INTRODUCTION 


The Naval Civil Engineering Laboratory (NCEL) has developed a line of semipermanent 
mooring and anchor devices called propellant-embedment anchor (PEA) systems. They come 
in a family of sizes of 10K, 20K, 100K, and 300K, where K denotes 1,000 pounds of nominal 
pullout capacity in conventional seafloor soil. The assembled PEA system is composed of a 
launcher and a projectile. To deploy the system it is lowered over the side of the ship, and when 
the ocean floor is sensed, the propellant charge is detonated, and the launcher ballistically drives 
the projectile into the ocean floor. Figure 1 shows the sequence of events for conventional 
applications. The launcher assembly is composed of a modified gun barrel and a hemispherical 
reaction vessel. The reaction vessel provides recoil resistance to the downward thrust of the 
projectile caused by the propellant bum inside the modified gun barrel. The projectile is 
composed of a piston and an anchor fluke. The anchor flukes come in varied configurations to 
facilitate penetration into different ocean floor materials -- mud, clay, silt, sand, gravel, coral, 
or rock. In this study, the material of interest is rock. 



Figure 1 

Deployment of a propellant-embedment anchor system into conventional seafloor soil. 
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Penetration of rock anchor projectiles depends on parameters such as initial kinetic 
energy, geometry of the projectile, interface friction, and constitutive properties of the seafloor 
rock material. Given these parameters, it should be possible to construct simplified and reliable 
models that predict the depth of penetration and the residual stress distribution surrounding the 
embedded anchor. The pullout resistance may be predicted once the confinement stress and the 
interface friction are known. 

Data from various analytical and numerical investigations of the problem need to be 
developed to form the basis for simplified penetration and extraction model development. These 
simplified models, employing physical as well as empirical and interpolated parameters, would 
be used to predict the performance of PEA systems in seafloor rock media, including pullout 
capacity. 

Based upon full-scale field tests in ocean sandstone and basalt rock, Beard (1984) 
concluded that the technology of PEA is feasible. Penetration depths up to one meter and 
maximum extraction forces of 200,000 pounds were obtained using the 20K propellant 
embedment anchor system. However, the penetration depths could not be correlated with anchor 
holding capacities, as concluded by Wadsworth and Beard (1980). 

Therefore, the development of analytical models and numerical methods to simulate 
anchor penetration and extraction processes was initiated to improve understanding of the 
penetration and extraction processes. This also includes studying the response of seafloor rock 
materials under high rates of loading, and developing validated confinement and friction models 
for describing the anchor extraction process, both of which are being pursued elsewhere. 

Analytical and numerical solution approaches to penetration problems have concentrated 
on predicting the depth of penetration and have, for the most part, ignored the structural response 
of the penetrator. Most approaches assume the penetration process can be formulated as 
axisymmetric with a rigid projectile. Deformations of a steel anchor projectile are often small 
compared to the target material, but stresses in the projectile may be very large. 

Beard (1984) recorded some mechanical failures of the PEA systems, such as kinking and 
birdcaging of the mooring cable as well as shackle fractures, in the field tests. These failures 
were attributed to high deceleration forces. Shugar and Wager (1982) described structural 
failures of the larger PEA reaction vessel due to overpressure from the gun barrel muzzle. 
Therefore, in addition to prediction of the penetration process, prediction of the structural 
dynamics of the PEA system, including the projectile, is necessary since it imposes structural and 
mechanical constraints on system performance. Further, structural dynamics calculations can 
assist in the evaluation of alternative shapes of projectiles which are subject to mechanical 
constraints. Consequently, when the survivability of the anchor projectile is important, 
deformable projectile models are required. For these reasons, a deformable projectile model is 
demonstrated in this study instead of a simple rigid projectile model. 

Objective 

The objective of this preliminary study was to develop a numerical modeling technique 
for understanding the residual stresses surrounding deformable anchor projectiles that have been 
ballistically embedded in rock media. The approach was to include explicit finite element 
technology to model the nonlinear dynamic response of the projectile and surrounding rock 
media. 
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Further, the effects of various numerical model parameters on depth of penetration and 
confinement stress were to be studied to promote development of hypotheses for simplified 
models of the penetration process. 

Scope 


The explicit finite element computer program DYNA3D (Whirley and Hallquist, 1991) 
was employed to model the nonlinear structural dynamics of the penetration problem. The three- 
dimensional dynamic finite element solution approach included nonlinear material behavior and 
slide line interface behavior. 

A static spherical cavity expansion theory for brittle materials (Forrestal and Longcope, 
1990) was also pursued for general background into analytical procedures used for predicting 
penetration processes. Another potentially more powerful analytical solution was also pursued 
in the hope that it could be used to independently verify at least some part of the numerical 
solution such as prediction of penetration depth. A dynamic cylindrical cavity expansion theory, 
presented by Zhou (1988), was investigated in this regard. 

The publications of Schwer and co-workers, Schwer, et al. (1988) and Schwer and Day 
(1991), were very useful, and the first paper was influential in the direction of the present 
research. It discussed the Lagrangian and Eulerian numerical frameworks for penetration 
problems and concluded that the former framework is more advantageous, particularly for 
computing the structural response of the penetrator. Further, the use of a prescribed artificial 
tunnel was recommended in the construction of the numerical model of the target medium to gain 
control of simulation time at the expense of having to assume a predetermined penetration path. 
The first publication also addresses axisymmetric analysis with DYNA2D while the second 
publication addresses three-dimensional analysis with DYNA3D. The latter publication became 
known to us only after the present study was completed. It presents three-dimensional studies 
of penetration and perforation processes in concrete and steel slabs. Where a single plane of 
symmetry exists in such problems, the use of a technique called the "local modified symmetry 
constraint" method, in lieu of the tunnel technique, is recommended. However, the tunnel 
technique, which is advocated in the present study, is currently the only alternative when fully 
three-dimensional simulation is required. For example, when the effects of in-situ rock fracture 
conditions are to be included, a fully three-dimensional model will be required to handle material 
anisotropy. 

Zhou (1988) and Goldsmith and Zhou (1989) have demonstrated some agreement between 
DYNA2D predictions and laboratory measurements of projectile penetration into rock specimens. 
In the present study, a measure of verification was also sought by comparing Zhou’s DYNA2D 
results with DYNA3D results for an axisymmetric model. 

Material properties encountered in a literature survey were used in this study to specify 
material behavior of the rock and projectile. Material properties for rock media were found to 
vary significantly from publication to publication. This variation defines parametric ranges over 
which future numerical penetration models may eventually be exercised to sort out the 
consequences of the variation. For purposes of this study, the mechanical property data compiled 
by Goldsmith and Zhou (1989) were used in conjunction with elasto-plastic models of material 
behavior in the development of the numerical penetration model. 

The DYNA3D finite element program and the corresponding pre- and post-processor 
codes, INGRID and TAURUS, were installed on a Sun 4/260 computer workstation. The 
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computer simulation times reported are based on the processor speed for this hardware; the 
processor ran at 25 MHz and is rated at 11 mips or 4 MFLOPS. 


REVIEW OF PREVIOUS SOLUTION METHODS 
Analytical Methods 

Analytical methods for solving projectile penetration problems are based on either the 
cavity expansion theory or the differential area force law as explained by Heuze (1989). These 
methods seek closed lorm or analytical solutions for the pressure on the projectile tip and the 
depth of penetration. 

Cavity expansion theory employs either a spherical or a cylindrical geometry, "'"’he former 
is suited to penetration by blunt nose projectiles, while the latter is suited to penetration by 
conical nose projectiles. These theories may further be categorized as static or dynamic. 
Resistance to penetration is assumed to be due to the shear strength of the target material. 
Closed form solutions for pressure on the projectile tip and depth of penetration are provided. 
However, they do not predict the confinement stress on the projectile body which is needed for 
determining anchor extraction force. Further, they do not address the stress response of the 
projectile. 

The differential area force law provides explicit formulations for the normal stress and 
tangential stress at every point on the external surface of the penetrator. This approach is three- 
dimensional and consists of a six-degree-of-freedom analysis of a rigid body penetrator. The 
method is primarily empirical because the material coefficients are empirically determined. 
According to Schwer et al. (1988), this theory is the basis for the PENC02D program 
(Creighton, 1982), and deceleration histories of rigid body peneTators are easily and adequately 
predicted by this program for normal impacts. Clearly, PENC02D was developed to address 
target penetration, excluding the structural response of the penetrator. 

The application of static spherical cavity expansion theory for determining resistance to 
penetration into brittle materials, including cracking, has been pursued by Forrestal and 
Longcope (1990). This approach to an otherwise complex problem is compelling in its 
simplicity. It is useful since it can provide the stress acting on the tip of the projectile in the 
presence of elastic, plastic, and cracking material behavior in the target. In an auxiliary 
penetration model, the dynamics of a rigid spherical penetrator is considered, and the stress 
found by application of a spherical cavity expansion theory is integrated over the surface of the 
projectile nose to get the force resisting penetration. The analytical model helps in defining the 
parameters that are important to the determination of target resistance. Moreover, it serves as 
a guide to interpretation of results from numerical models as is demonstrated . ‘h; section 
"Analysis of Confinement Stress" in this report. 

Appendix A contains a review and derivation of the associated equations of the spherical 
cavity expansion theory developed by Forrestal and Txmgcope (1990). For brittle materials the 
shear strength varies with confinement pressure. Thus, the radial stress which will retard 
penetration of a hemispherical projectile will vary with the pressure dependent shear strength 
parameter for the target material. This dependency is shown in Figure 2 for both a ceramic and 
plain concrete target material. The increase in hardness of the target due to the assumption of 
incompressible material behavior is also indicated for each material. This, for example, suggests 
the possibility that saturated soil or rock materials on the seafloor would develop a similar 
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Figure 2 

Radial stress from spherical cavity expansion theory. 


increase in apparent hardness since water in the pores would not have time to migrate away from 
the projectile during the highly transient event. 

Zhou (1988) presents a full discussion of a dynamic cylindrical cavity expansion theory 
for anchor penetration problems. Though simplifications designed to facilitate application were 
introduced, Zhou apparently felt the theory still remained too cumbersome for application and 
did not apply it in conjunction with his numerical studies. The set of equations derived were not 
exercised by Zhou, but they potentially offered a means in the present study to independently 
validate numerical predictions of depth of penetration using DYNA3D. This was an inducement 
to further analysis of these equations in the present study. However, it was soon found that the 
expression given for shear resistance contained an infinite series which diverged. The 
corresponding equations and the nth-term test for identifying divergent series are presented in 
Appendix B. The cause of the errant equation is unknown, and further work is required to apply 
this theory. 

DYNA2D Studies 

Zhou (1988) investigated the normal impact penetration of a fastener projectile into 
concrete and sandstone. He used theoretical, experimental, and numerical approaches in the 
study. Comparisons obtained between numerical and experimental data were reasonably good. 
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As mentioned, no comparisons were made using the cylindrical cavity expansion theory approach 
though an extensive presentation of the theory was given. 

Zhou and Goldsmith (1989) also investigated factors related to PEA projectiles embedded 
in rock. They attempted to provide an analytical solution which was directly applicable to PEA 
projectiles based on the theory of elasticity, but certain nonlinear terms defied solution and the 
analytical effort was terminated. They also reported on a preliminary numerical study. The 
results from this numerical study were compared to experimental results with reasonably good 
correlation. Unfortunately, the DYNA2D input data file necessary for repeating the analysis was 
not presented and replication of the analysis using DYNA3D, as described herein, was therefore 
not necessarily precise. 

Numerical modeling techniques based upon using DYNA2D to describe the target were 
emphasized in Schwer, et al. (1988). Success rested on the prescription of a small tunnel at the 
center of the axisymmetric model of the target which also predetermines the path the penetrator 
would be constrained to follow. The tunnel had the practical modeling effect of reducing the 
severity of element distortion in the finite element mesh, and reducing the computation time 
without apparently otherwise affecting the results. The agreement with experimental results was 
considered acceptable. 


THREE-DIMENSIONAL FINITE ELEMENT MODELS FOR PENETRATION 

Numerical modeling of the penetration process can be approached by using finite element, 
finite difference, or discrete element methods. Both Lagrangian and Eulerian formulations of 
finite element and finite difference methods have been used, and each has its advantages and 
limitations. 

Eulerian codes are generally used for fluid dynamic problems where the fluid enters and 
exits from a control volume. They are also applicable to large strain problems in solid 
mechanics. They have been applied to the penetration problem with the target being the control 
volume. They are more advantageous than the Lagrangian codes when very large strains are 
involved because no distortions can occur to the control volume grid; this obviates the need to 
re-zone the model grid at various times during the simulation of the problem. The difficulty with 
using Eulerian codes involves inaccurate descriptions and calculations of material boundaries and 
interface conditions. Processing structural response information from Eulerian codes requires 
careful treatment of data on how the penetrator is flowing through the grid (Schwer et al., 1988). 
The assumption of a rigid penetrator often accompanies an Eulerian code application to remove 
this modeling difficulty. If a rigid penetrator model is satisfactory, Eulerian codes may be 
competitive for penetration problems. 

Lagrangian codes require the model grid to follow displacements of the material. They 
provide a more accurate prediction of interface definition and behavior as well as penetrator 
structural response, if deformations of the grid can be managed so that re-meshing is avoided or 
minimized. The three-dimensional program DYNA3D is strictly a Lagrangian finite element 
code, while the two-dimensional program DYNA2D is primarily a Lagrangian code with limited 
Eulerian formulation capability. 

Understanding how the anchor extraction force is dependent on the damaged rock material 
adjacent to the embedded projectile is crucial in predicting performance of anchor projectile 
systems as noted by Wilson et al. (1990). Predicting the resulting post anisotropic state of the 
target material, as well as including the in-situ anisotropic state of the target, are potentially 
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important features of successful numerical simulation. As these features are not generally 
axisymmetric, three-dimensional models are required to replicate them. Therefore, application 
of DYNA3D vice DYNA2D for simulation of the penetration process has been pursued in the 
present study to promote procedures and techniques for calculating fully three-dimensional 
confinement stresses surrounding embedded projectile anchors. DYNA3D is also currently 
regarded as being better supported for material models that can be used for describing the 
behavior of rock. Additional advantages of DYNA3D vice DYNA2D for anchor projectile 
penetration simulation might also include the potential for any future extension to oblique angles 
of impact and non-rotationally symmetric anchor projectiles. 

Three-dimensional modeling in the present investigation employed quarter models of 
axisymmetric configurations because they are sufficient to develop three-dimensional modeling 
techniques while minimizing simulation time. Fully three-dimensional model demonstrations 
were beyond the scope of this initial investigation. 

The critical feature of the DYNA codes relative to the objectives of the present project 
is their capability for handling sliding interfaces (Hallquist, 1978). The numerical treatment 
involves slave and master node lines defined by the user such that a slave node is constrained to 
slide on a line between two adjacent master nodes subsequent to impact until a tensile force 
develops. A separation between a slave and master line is referred to as a void in DYNA3D 
parlance. The designation of the master surface is apparently arbitrary (Hallquist, 1977). The 
dynamics of contact/impact and masicr/slave interfaces as well as other pertinent features, are 
discussed further in Appendix C, where a brief description of DYNA3D is provided. 

Previous DYNA3D Studies 

Rainsberger Example. Rainsberger (1988) developed an example of a three-dimensional 
finite element penetration model using DYNA3D. He fully described the necessary parameters 
needed to reproduce the computational results for a projectile penetrating a rock target. The grid 
technique for both the target and the projectile was described. 

The Rainsberger example was reproduced in this study at the outset to verify that the 
software system consisting of INGRID, DYNA3D, and TAURUS was operating as intended. 
INGRID and TAURUS are, respectively, pre- and post-processor programs for DYNA3D. The 
example consisted of a two-part steel projectile penetrating a target. The target material behavior 
was modeled by a soil/crushable foam material model (DYNA3D Material Type 5). All the 
parameters were defined and presented in an INGRID input data file. The example contained 
a quarter model (two planes of symmetry were exploited) with three slide surfaces defined as 
shown in Figure 3. Slide surface 1 was used to tie the two-part steel projectile together. Slide 
surface 2 was a sliding interface (including separation or void formation) between the projectile 
and the target. Slide surface 3 was used to tie the differently zoned target model together and, 
thus, avoid the requirement for special transition elements at the interface. After the skill 
necessary for using the pre- and post-processor programs and the slide line capability was 
developed, the example was reproduced as published. Though not shown here, the resulting 
pressures in the target, on the projectile, and at the projectile/target interface were typical of the 
type of data sought in the present project. The Rainsberger example thus provided the initial 
guidance required to apply DYNA3D to the problem at hand. 
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Figure 3 

Rainsberger penetration model and location of slide surfaces. 


Schwer and Day Simulations. Schwer and Day (1991) have employed DYNA3D for 
oblique angles of penetrator impact into concrete targets where the penetration path is 
predetermined using the tunnel method. They also recommend a method called the local 
modified symmetry constraint technique, where penetration may take any path restricted to a 
predetermined plane of symmetry within an otherwise arbitrary three-dimensional target medium. 
(Both the target geometry and target material behavior must be symmetric with the plane of 
symmetry.) The method proports to work similarly to the tunnel method, in effect exploiting 
DYNA3D features for releasing symmetry constraints and providing slide surfaces. The 
simulations presented are remarkable and fully demonstrate the capability of DYNA3D for this 
class of problem. 

Fastener Model Using DYNA3D 

A published DYNA2D analysis of fastener projectile penetration by Zhou (1988) was 
repeated using DYNA3D. The small projectile, shown in Figure 4, was 4 millimeters in 
diameter, 35 millimeters long, and had a 2.5 caliber radius head with an ogival nose. 
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Figure 4 

Size of small fastener projectile in millimeters. 

The three-dimensional version of Zhou’s axisymmetric model contained a frictional sliding 
interface with voids or gaps permitted to form between the projectile and target. A quarter 
finite element model was constructed exploiting two planes of symmetry of the axisymmetric 
configuration. Correspondingly, two planes of symmetry were used in the INGRID input data 
file to effect the quarter model. The typical INGRID input data file generated for this simulation 
is presented in Appendix D. The file has been annotated for convenience to the reader as to the 
meaning of the INGRID commands used. Additional data required for the projectile penetration 
model had to be reconstructed from Zhou’s narrative description since actual input data files for 
his DYNA2D model were not included in his publication. Consequently the recreated DYNA3D 
model may not exactly replicate his model. 

The projectile was modeled with a small artificial bore hole (0.4 mm in diameter) to 
effect the formation of a hexahedron element when the plane containing the lines were rotated 
to make a solid. This is crucial to numerical efficiency of the overall modeling technique 
because it leads to control of the time step size. The artificial bore does not affect the results 
since its surface does not come in contact with the target tunnel. 

Figure 5 shows the two-dimensional projectile geometry used to generate the three- 
dimensional geometry for the model. The designated line numbers correspond to numbered line 
definitions in the INGRID input data file (see Appendix D). This two-dimensional geometry is 
then rotated about the axis of symmetry to form solid elements. The projectile was divided into 
a grid 2 elements wide by 24 elements long and then rotated 90 degrees in five successive 
increments or copies to create 240 eight-node solid elements. Figure 6 shows the initial two- 
dimensional geometry of the target and projectile together. 

Zhou’s target was a cylindrical block of rock having both a height and diameter of 150 
millimeters. In the corresponding model, the first 15 millimeters were finely zoned to 
accommodate the large amount of deformation expected in the region. The planar mesh in this 
zone was 4 elements wide by 102 elements deep and was rotated in five successive increments 
to form a quarter model of the inner part of the target. The second or outer part of the target 
was divided into a 12-element-wide mesh with the elements increasing in width with distance 
from the center. The planar mesh was 34 elements deep and was rotated in five successive 
increments to form a quarter model of the outside solid portion of the target. The two different 
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Figure 5 

The two-dimensional geometry of the projectile from Zhou’s penetration 
problem. (The line numbers refer to the line definitions in the 
INGRID input data file. The scales are in millimeters.) 

meshes of the target were joined using a tied sliding boundary. No sliding or separation occurs 
at this boundary; this use of the slide boundary capability in DYNA3D is a convenient and 
expedient way to avoid dealing with transitional elements in preparing the model. A total of 
4,080 eight-node solid elements make up the target model. 

The resulting three-dimensional finite element model of the target and projectile was 
composed of 4,380 eight-node solid elements, and is shown in Figure 7. The bottom nodes of 
the target model were restrained against displacement. 

A small artificial tunnel (1.4 mm diameter) was introduced at the centerline of the target. 
The tunnel reduced the risk of extreme distortion of the elements in the target. This in turn 
facilitates the use of the Lagrangian approach for describing material deformation by avoiding 
having to remesh the model. Once introduced, the size of the tunnel will significantly affect the 
computer time required; a large tunnel diameter uses less computer time than a small tunnel 
diameter. This is because explicit finite element computer programs such as DYNA3D cannot 
step forward in the time integration procedure for solving the equations of motion any faster than 
the smallest element allows. For the integration to remain stable, the time step size must remain 
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Figure 6 

Two-dimensional representation of the projectile and target from Zhou’s 
penetration problem. (The line numbers refer to the definitions 
in the INGRID input file. The scales are in millimeters.) 


smaller than the transit time of an elastic wave traversing the smallest element in the model. 
However, the effect of the size of the artificial tunnel on computational results is of concern and 
was also addressed by Schwer et al. (1988). It is important to determine this effect and it will 
be addressed later in this report as well. 

A sliding interface was introduced between the areas of contact of the target and 
projectile. This interface was assigned a frictional value of zero (the value used by Zhou), and 
gaps or voids were allowed to form during the simulation. 

Material properties for rock were sought from the literature to support development of 
the numerical model. Published values for rock media vary significantly. Some published 
mechanical properties for rock are compiled in Table 1. Soft and hard sandstone, and vesicular 
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and intrusive basalt materials are included. 

The data developed by Goldsmith and 
Zhou (1989) were used in the numerical 
penetration model in the present study. 

The material properties used for modeling 
the steel projectile anchors are listed in 
Table 2. 

The steel projectile material 
behavior and rock target material behavior 
were both modeled using DYNA3D 
Material Model #3, an elastic-plastic 
material model with isotropic or kinematic 
hardening. Specific parameter values 
used for the projectile and target material 
models are given in Tables 3 and 4, 
respectively. As nearly as can be 
determined, these values were used in the 
axisymmetric model by Zhou, who used 
DYNA2D’s version of Material Model 
#3. The uniaxial stress-strain curve 
shown in Figure 8 depicts elastic-plastic 
material behavior for kinematic hardening (0 = 0) and isotropic hardening (0 = 1.0). Kinematic 
hardening was assumed for both the projectile and the target. 

A typical set of data plots and graphs of results from the simulation is shown in Figures 
9 through 11. Appendix E contains information on a sample session for pre-processing, 
executing, and post-processing a projectile anchor penetration simulation using, respectively, 
INGRID, DYNA3D, and TAURUS. The TAURUS batch file used to create the plots and 
graphs in Figures 9 through 11 is given at the end of Appendix E. The INGRID input data file 
can be found in Appendix F where it is labeled as try#18. 

The results from the DYNA3D model were near the expected values for the maximum 
depth of penetration. Zhou reported a depth of penetration of 20 mm as shown in Figure 12(a). 
In both cases, the initial impact velocity was 100 m/s. The DYNA3D simulation computed a 
maximum depth of 15 mm as shown in Figure 10(a). This represents a 25 percent difference 
between the two simulation models. In a displacement-based finite element formulation, the 
exact displacement solution is always approached monitonically from below, though reduced 
integration (see Appendix C) mitigates this argument. In this case, the DYNA2D simulation 
must be regarded as more accurate since the configuration modeled is axisymmetric. Therefore, 
that the DYNA3D displacement is less than the DYNA2D displacement is somewhat as expected. 

Other differences between the DYNA2D and DYNA3D simulation results involved the 
amount of rebound and final velocity of the projectile. The reported velocity for the DYNA2D 
simulation reversed and approached a zero residual value, and the reported rebound of the 
projectile was less than 10 mm. The velocity for the DYNA3D simulation also reversed but 
approached a nonzero steady state residual value causing the projectile to back out of the cavity. 
This was apparently due to residual cavity pressure acting on the ogival nose of the projectile and 
the absence of friction. 



Figure 7 

DYNA3D finite element penetration model 
for fastener projectile. 
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Table 1 

Mechanical Properties of Rock Materials 


Reference 

Source 


(kg/m 3 ) 


Young’s 

Modulus 

(GPa) 


Beard (1984) 

2231 

Stowe (1984) 

2291 

Goldsmith & Zhou 
(1989) 

2164 

Zhou (1988) 

2130 

Brunette & Goldsmith 
(1989) 

1940 



Poisson’s 

Ratio 

Porosity 

(%) 

Ultimate 

Compressive 

Strength 

(MPa) 

Ultimate 

Tensile 

Strength 

(MPa) 

Soft Sandstone 

- 

- 

10.13 

- 

0.35 

20.7 

15.24 

0.161 

- 

- 

23.37 

- 

0.20 

18.0 

45.00 

4.55 

- 

- 

7.05 

- 


Rock Type: Hard Sandstone 


Wadsworth & Beard 
(1980) 

2600 

26.4 

0.01 

Goldsmith & Zhou 
(1989) 

2588 

47.34 

- 

Brunette & Goldsmith 
(1989) 

2500 

4.98 

- 



Rock Type: Vesicular Basalt 


Wadsworth & Beard 
(1980) 

2167 

7.29 

0.167 


30. 

33.13 

Beard (1984) 

2325 

26.24 

- 


- 

17.31 

Goldsmith & Zhou 
(1989) 

2248 

26.24 

- 


26. 

40.4 

Brunette & Goldsmith 
(1989) 

2970 

8.29 

- 


- 

123. 



Rock Type: 

Intrusive Basalt 



Beard (1984) 

2710 

- 

- 


- 

50.48 

Stowe (1984) 

2662 

34.48 

0.15 


7.0 

87.08 

Goldsmith & Zhou 
(1989) 

2798 

51.88 

- 


- 

99.93 

Brunette & Goldsmith 
(1989) 

2500 

3.53 

- 


- 

56.3 




































Table 2 

Material Properties of AISI 4340 Steel 


Density . 

. 7832 kg/m 3 

Young’s Modulus .... 

. 200 GPa 

Yield Strength . 

. 310-793 MPa 

Tensile Strength . 

. 552-2000 MPa 


Table 3 

Projectile Material Properties for Material Model #3 


Density 

7.8E-6 kg/mm 3 

(7800 kg/m 3 ) 

Young’s modulus, E 

209E6 kg/(mm-s 2 ) 

(209 GPa) 

Poisson’s ratio 

0.3 


Uniaxial yield stress, a Q 

1.2E6 kg/(mms 2 ) 

(1.2 GPa) 

Tangent modulus, Ej 

1E6 kg/(mms 2 ) 

(1.0 GPa) 

Hardening parameter, 0 

0 



Table 4 

Target Material Properties for Material Model #3 


Density 

2.13E-6 kg/mm 3 

(2130 kg/m 3 ) 

Young’s modulus, E 

13.62E6 kg/(mm s 2 ) 

(13.62 GPa) 

Poisson’s ratio 

0.2 


Uniaxial yield stress, a D 

0.026 kg/(mm s 2 ) 

(0.02 GPa) 

Tangent modulus, Ej 

0 


Hardening parameter, 0 

0 
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= maximum stress 
E = Young’s modulus 

o o = yield stress 

E t = tangent modulus 

P = hardening parameter 

e = strain 

a = effective stress 


Figure 8 

Uniaxial stress-strain curve showing elastic-plastic material 
behavior for kinematic hardening 03=0) and isotropic hardening 03=1) 
(Hallquist and Whirley, 1991). 


Without friction between the projectile and target, and with a zero value prescribed for 
the target hardening modulus, the problem as modeled is physically unrealistic. The noted 
differences in results may also be due to misinterpretation of Zhou’s model. Again, the results 
of the DYNA3D model should be expected to approach those for the DYNA2D model as the 
three-dimensional mesh is refined around the perimeter of the model. In this respect, the 
DYNA3D model must be viewed only as an approximation to the DYNA2D model. The 
differences noted between the results of these two simulations otherwise are not explained. 

The computer processing time for a 0.5-ms DYNA3D simulation on the Sun 4/260 
workstation was approximately 45 hours. The amount of computer time on this simulation is 
excessive for two reasons. First, the projectile had reached a steady state condition within 
approximately 0.25-ms. Thus, the simulation could have been terminated with no loss of useful 
data much earlier. Second, the number of elements in the target model was excessive. This was 
because most of the elements that experienced significant pressure were within 5 projectile radii 
of the centerline (Figure 11). A similar result was noted by Schwer et al. (1988). The finite 
element model needlessly extended well beyond this range. It is therefore surmised that a 
satisfactory simulation could have been obtained by halving the duration and using a substantially 
reduced finite element model of the target. 
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(b) Side profile. 


Figure 9 

Projectile embedded in the target at end of the simulation. 
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Figure 11 

Contours of residual pressure in the tarj 
(Units for the pressure are kg/(mm*s 


(a) Isometric view. 


(b) Top view. 


(c) Side view. 
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Axial Velocity (cm/ps) Penetration Depth (cm) 



(a) Displacement. 



(b) Velocity. 


Figure 12 

Projectile response from DYNA2D (Zhou, 1988) 
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NUMERICAL MODEL TECHNIQUE AND DEMONSTRATION 


The basic DYNA3D numerical model discussed in the previous section was exercised 
further to evaluate the effects of various modeling parameters. The parameter study would help 
verify the numerical model by investigating areas of uncertainty that were raised during the 
preliminary analysis and the literature review, such as artificial tunnel size, effect of friction, 
material model selection, and amount of strain hardening. After the parameter study, the 
geometry for the model was changed to reflect a larger anchor projectile for rock anchor 
penetration simulation. 

Parameter Study 

Effect of Target Tunnel Diameter. The effect of artificial tunnel size was addressed 
by Schwer (1988) in DYNA2D models. An upper limit on the tunnel diameter which minimizes 
computation time and maintains accuracy is perhaps even more meaningful for DYNA3D 
applications than for DYNA2D applications because of higher compute times in three 
dimensions. A suitable range for the tunnel diameter was sought which would not compromise 
simulation results and would minimize computation time. The diameter of the tunnel was varied 
in five runs to determine how it affects the computed response. (In Appendix F these runs 
correspond to INGRID input data files labeled try #7, 8, 9, 15, and 16.) The hardening modulus 
used for these simulations was one percent of Young’s modulus, and the coefficient of friction 
used between the projectile and the target was 0.1 as suggested by Goldsmith and Zhou (1989). 

Figure 13 shows the computed response of the projectile for various tunnel diameters. 
The displacement histories show various amounts of projectile rebound associated with the 
various tunnel diameters. The smallest diameter causes the largest rebound of the projectile. 
Rebound was not discussed in the previous DYNA studies reviewed. However, the amount of 
rebound is of interest because the residual penetration depth can clearly play a significant role 
in predicted anchor holding capacity. The projectile velocity neared a steady state velocity of 
20 m/s upward for the smallest tunnel case, and was still rebounding when the simulation was 
manually terminated. For all the larger tunnels, the projectile came to rest at various rebound 
depths. 

The acceleration history also varied significantly with target tunnel diameter. Mechanical 
failures of anchor systems have been associated with high deceleration rates. Thus, these data 
are important to development of modeling technique. Larger diameter tunnels result in lower 
computed values for projectile deceleration. Also, since the 0.40-mm and 1.02-mm tunnel 
simulations are similar, a lower limit for the tunnel diameter in this range is suggested at which 
higher deceleration will not occur for the specific problem modeled. 

The contours of pressure shown in Figure 14 provide insight into the fundamental cause 
of projectile rebound. These data are shown at 0.5 ms into the simulation, the point where the 
projectile generally had come to rest. They may be regarded as residual pressures in the target 
material. The largest positive (compressive) pressures occur near the ogival nose of the 
projectile, where force is mobilized to push the projectile back out of the cavity. Smaller 
diameter tunnels lead to larger contact areas on the nose of the projectile over which the residual 
pressures act, and, thus, a larger force is mobilized to push the projectile back out of the cavity. 
The deformation of the tunnel is not as pronounced in Figure 14(e) and the force mobilized is 
small when compared to the case shown in Figure 14(a) where a well-pronounced cavity has 
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Figure 13 

Effect of target tunnel diameter on penetration depth, velocity, and acceleration. 
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Figure 14 

Residual pressure contours in the target for various tunnel diameters. 
(Units for contour values are kg/(mm*s 2 ) 
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formed. The force causing rebound is eventually offset by the friction force which opposes 
rebound, and the projectile finally comes to rest in the case of the three larger tunnels. 

The effect of the tunnel diameter on penetration depth was non-dimensionalized by 
dividing the tunnel diameter by the projectile diameter to facilitate a useful definition of an upper 
limit for the tunnel diameter. Figure 15 shows the effect of the ratio on the maximum 
penetration depth. As the ratio approached zero, the slope of the curve approached zero. These 
data indicated that a reasonable upper limit for the ratio was 0.25. Consequently, as long as the 
artificial tunnel diameter remains one quarter or less of the projectile diameter, the DYNA3D 
penetration model should be sufficiently insensitive to the presence of an artificial tunnel. The 
ratio of 0.25 appears to be an optimal balance between accuracy and computation time for 
DYNA3D simulations of anchor penetration into rock materials. Schwer and Day (1991) report 
that a ratio of 0.4 or less does not affect the penetrability of the target model when using 
DYNA2D. That standard applied to the present DYNA3D simulation would lead to a 20 percent 
error in predicted penetration depth. Thus, the critical target-to-projectile diameter would seem 
to depend on whether DYNA2D or DYNA3D is being used. 


20.0 

18.0 

16.0 

14.0 

12.0 

10.0 

8.0 

6.0 

4.0 

2.0 

Penetration Depth (mm) 





































. 






( 

>— - 

—e-—' 






































u.« 

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 

Target Tunnel Diameter/Projectile Diameter 


Figure 15 

Normalized chart for the effect of the target/projectile diameter 
ratio on the maximum penetration depth. 
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Effect of Friction Coefficient. The effect of friction at the projectile/rock interface was 
studied to develop skill in applying the friction slideline capability in DYNA3D, and to 
determine the sensitivity of the model to the friction coefficient. Since little is known about the 
coefficient of friction for rock anchor penetration, and coefficient of friction values vary 
significantly in the literature, the model’s response to different values needed to be quantified. 
The effect of the friction coefficient prescribed for the sliding interface between the projectile 
and target was varied from zero to 0.1 in three simulations. The tunnel-to-projectile diameter 
ratio was held constant. The resulting effect on penetration model response is shown in 
Figure 16. Predictably, a smaller coefficient of friction results in a larger penetration depth. 
The variation in maximum penetration depth was approximately 20 percent. However, this 
diminished to 10 percent or less when rebound was considered and residual penetration depth was 
the basis of comparison. 

The residual rebound velocity for the frictionless case and for the 0.05 coefficient of 
friction case was nonzero. Had the simulation been allowed to continue, the projectile would 
have exited the cavity for the frictionless case and perhaps as well for the 0.05 case. These data 
indicate that only simulations with a coefficient of friction equal to 0.1 or greater would prove 
useful for anchor extraction studies because static residual conditions are essential to such studies. 
The residual pressure contour data for this simulation are shown in Figure 14(c). 

The data also show that larger values of interface friction result in higher maximum 
deceleration levels as expected. Further, greater values of interface friction show lesser 
oscillation in the projectile deceleration graphs. This is perhaps because interface friction reduces 
relative (tangential) motion between the projectile and rock, and the Type 3 sliding surface 
algorithm is invoked less often. In doing so, the response of the projectile is smoother. It would 
be expected that the frequency content of the impact time histories would be higher for the higher 
friction simulations shown because these impacts are sharper. 

Effect of Different Material Models. Various DYNA3D material models were 
considered in search of a model that would suitably describe material behavior of the rock target. 
A good material model should allow for the cracking tendency of rock. DYNA3D contains four 
potential material models for modeling rock targets: kinematic/isotropic elastic-plastic model 
(#3), a soil and crushable foam model (#5), a concrete/geological model (# 16), and a geologic 
cap model (#25). 

In the previous parameter studies, DYNA3D material model #3 was used exclusively. 
Here the rock target material model was changed to a soil and crushable foam model, DYNA3D 
material model #5, which can also be used to model rock including tensile cracking. The 
parameter values used for this model are given in Table 5. These data prescribe a pressure 
versus volumetric strain graph of the form shown in Figure 17. The values for the yield function 
constants were obtained from the DYNA2D study by Zhou (1988) which intended to replicate 
sandstone targets. Material model #3 was again used for the projectile (Table 3 and Figure 8). 

Unfortunately material model #5 failed to provide a useful penetration model with a 
sandstone target when using the parameter values described. As shown in Figure 18, the target 
material displaced away from the projectile. These results would not be useful in conjunction 
with anchor extraction studies because no condition of confinement was predicted. A comparison 
of the penetration model response for material models #3 and #5 is shown in Figure 19. Overall, 
the results are similar. Material model #5 constituted a softer target and allowed approximately 
10 percent more penetration, and the projectile reached a small positive steady state velocity and 
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Figure 16 

Effect of friction coefficient on projectile penetration depth, 
velocity, and acceleration. 
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Figure 17 

Pressure versus volumetric strain curve for soil and crushable foam 
material model (Whirley and Hallquist, 1991). 


Table 5 

Target Material Properties for Material Model #5 


Density 

2.13E-6 kg/mm 3 

(2130 kg/m 3 ) 

Shear modulus 

5.68E6 kg/(mms 2 ) 

(5.68 GPa) 

Bulk unloading modulus, K u 

11.25E6 kg/(mms 2 ) 

(11.25E6 GPa) 

Yield function constant 0 

1.33 kg/(mms 2 ) 


Yield function constant 1 

60,270 kg/(mm s 2 ) 


Yield function constant 2 

0 


Pressure cutoff for tensile 
fracture, p cut 

0 


No. of points on the p vs. e v curve 

6 


Volumetric strain values, « v 

0 .0025 .01 .04 

.05 .055 

Corresponding pressure values, p 

0 15000 21250 34000 

41300 50000 










Figure 18 

Projectile embedded in target characterized by material model #5. 


was backing out of the cavity. The projectile acceleration response is very similar for material 
models #3 and #5. 

In further attempts with material #5, an error was encountered when trying to assign any 
value other than zero to the pressure cutoff for tensile fracture; a small negative value of ten and 
a small positive value of ten encountered run-time errors. DYNA3D material models #16 and 
#25 were not considered in this study. Future rock anchor penetration studies should investigate 
these models as well. 

Effect of Hardening Modulus. Using material model #3 and DYNA2D, Zhou (1988) 
reported a 25 percent reduction in penetration depth when using zero hardening instead of 0.5 
percent of Young’s modulus for the hardening modulus value. Therefore, the hardening modulus 
in material model #3 was varied to see if a similar reduction occurred in a DYNA3D simulation. 
Goldsmith and Zhou (1989) suggested using a nominal value of 1.0 percent of Young’s modulus 
for the hardening parameter. The hardening modulus (or tangent modulus, Ej, in Figure 8) 
clearly increases the resistance of the target, and a reduction in penetration depth is anticipated 
with higher values. It is expected to be a critical material model parameter in modeling 
penetration. 

The tangent modulus for this study was set at zero and 1.0 percent of Young’s modulus 
for two successive simulations. The material properties used for the target (see Table 4) were 
otherwise unchanged. The coefficient of friction between the projectile and target was 0.1. 
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Figure 19 

Effect of rock target material models on projectile penetration 
depth, velocity, and acceleration. 
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Figure 20 compares the penetration response for different values of target hardness. The 
penetration depth was about 25 percent less for the hardening modulus value of 1.0 percent of 
Young’s modulus compared to no hardening. The maximum deceleration of the projectile was 
correspondingly greater since the penetration depth decreased. Corresponding residual pressure 
contours in the rock target are shown in Figure 21. These data show that the average direct 
residual stress (residual pressure) on the projectile is higher for the lower hardening modulus. 

Larger Anchor Projectile Study 

The penetration model for a larger anchor projectile was constructed from knowledge 
gained in replicating Zhou’s small fastener penetration model and from the aforementioned 
parameter study. The larger projectile was 2.54 cm in diameter and 17.78 cm long. The ogival 
nose is defined by the ratio between the nose arc’s radius and the projectile diameter; in this 
case, the caliber radius head was 3.0. A sketch of the projectile is shown in Figure 22. The 
following description of the construction of the DYNA3D penetration model for the larger 
projectile is included for completeness. 

Construction of the finite element model of the projectile was facilitated by including a 
small artificial bore 0.508 cm in diameter at the center. This prevented very small hexahedron 
solid elements from being formed in the development of the projectile model. The projectile was 
divided into a grid 2 elements wide by 14 elements long and then rotated 90 degrees about the 
projectile centerline in five successive increments to create 140 eight-node solid elements which 
constitutes a one-quarter solid model of the projectile. The material behavior of the steel 
projectile was modeled using an elastic-plastic material model (DYNA3D model #3). The 
parameter values used for this material model are given in Table 6. 

The target was modeled as a cylindrical block having a depth of 63.5 cm and a radius 
of 25.4 cm. To promote reasonable computer run times, a small artificial tunnel 0.635 cm in 
diameter was introduced along the centerline of the target. This corresponds to a tunnel-to- 
projectile diameter ratio of 0.25, which was found to be an optimum value in the preliminary 
study. The inner part (the first 10.16 cm along the radius) was finely zoned to accommodate 
substantial deformation gradients. The grid for this zone is 4 elements wide by 100 elements 
deep and was rotated in five successive increments to form a one-quarter solid model of the 
target. The outer part of the target was divided into a 3 by 25 mesh with the element size 
increasing by a factor of two for each increment from the center outward. The two different 
meshes of the target were joined using a tied sliding boundary. No sliding or separation was 
permitted at this interface. A total of 2,375 eight-node solid elements made up the one-quarter 
solid model of the target. 

The base of the target model was restrained against displacement. Target material 
behavior was modeled as an elastic plastic material (DYNA3D model #3) using the parameter 
values listed in Table 7. A nominal value of 1 percent of Young’s modulus was used for the 
tangent modulus as suggested by Goldsmith and Zhou (1989). 

A sliding interface was employed between the target and projectile. This interface was 
assigned a coefficient of friction value of 0.1 and permits development of gaps or voids during 
the simulation. The DYNA3D penetration model for the larger anchor projectile is constituted 
of 2,515 eight-node solid elements, and is shown in Figure 23. 
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Figure 20 

Effect of target hardening modulus on projectile penetration 
depth, velocity, and acceleration. 
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(a) 1 % of Young's modulus. 



Figure 21 

Contours of residual pressure in the target for varying hardening modulus values. 

(Units are kg/(mm*s 2 ) 
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Figure 22 

Larger anchor projectile. (Dimensions are in centimeters) 


Table 6 

Projectile Material Properties for 
the Larger Anchor Projectile 


Density 

7.80 g/cm 3 

(7800 kg/m 3 ) 

Young’s modulus 

2.00E12 g/(cm-s 2 ) 

(200 GPa) 

Poisson’s ratio 

0.3 


Uniaxial yield stress 

17.9E9 g/(cms 2 ) 

(1.79 GPa) 

Tangent modulus 

10.0E9 g/(cm s 2 ) 

(1.0 GPa) 
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Computer run times on the Sun 4/260 workstation rated at 11 mips were approximately 
18 hours, a reduction of over 50 percent in compute time from the DYNA3D model (Figure 7) 
of Zhou’s penetration problem. The reduction was achieved by halting the simulation when the 
projectile came to rest, and also by having reduced the number of elements in the model in the 
radial direction. Further reduction could be achieved by truncating the mesh in the vertical 
direction as well. 

In a subsequent simulation, the target material characterization was modeled with 
DYNA3D material model #5 so that the differences between the two material models for the 
target could again be assessed. The parameter values used for this model are listed in Table 8. 


Table 8 

Target Material Properties Used With Material 
Model #5 for the Larger Anchor Projectile 


Density 

2.16 g/cm 3 

(216 kg/m 3 ) 

Shear modulus 

8.44E10 g/(cm-s 2 ) 

(8.44 GPa) 

Bulk unloading modulus 

11.25E10 g/(cm-s 2 ) 

(11.25 GPa) 

Yield function constant 0 

.00133E-7 g/(cm-s 2 ) 


Yield function constant 1 

60.27 g/(cm-s 2 ) 


Yield function constant 2 

0 


Pressure cutoff for tensile 



fracture, p cut 

0 


No. of points on the p vs. € v curve 

6 

% 

Volumetric strain values, e v 

0 .0025 .01 

.04 .05 .055 

Corresponding pressure values, p 

0 15E7 21.25 

E734E7 41.3E7 50E7 


The projectile dynamic responses for material models #3 and #5 are shown in Figure 24. 
The initial impact velocity was 100 m/s. In this case, the projectile penetration was less for 
material model #5, which was opposite to the outcome for the smaller projectile (Figure 19). 
Otherwise the responses were similar. Moreover, material model ft 5 again failed to provide a 
residual stress field as shown by the absence of stress contours in Figure 25, due again to the 
target material having pulled away from the projectile. 
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Figure 24 

Effect of target material models on larger anchor projectile 
penetration depth, velocity, and acceleration. 
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(b) Material model #5. 


Figure 25 

Contours of target pressure for larger anchor projectile. 
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ANALYSIS OF CONFINEMENT STRESS 


The confinement stress field resulting from an embedded anchor projectile is prerequisite 
information for calculating extraction load-deformation behavior and the holding capacity of the 
anchor. The purpose of this section is to analyze the computed three-dimensional residual stress 
field that surrounds an anchor projectile, particularly at the interface where the stresses impinge 
on the projectile. Confinement stresses are first defined and discussed. Then computed results 
are discussed which demonstrate how to interpret confinement stress data in the context of anchor 
extraction studies. The effects of various parameters on confinement stress are also discussed. 

Description and Discussion of Residual and Confinement Stress Data 

In this section a distinction between residual and confinement stress is made. Residual 
stress is defined as the stress response computed at the time the anchor projectile has come to 
rest. Confinement stress is a subset of these data referred to the target elements lining the cavity 
formed by the impact. Specifically, the direct or normal residual stress component acting on the 
anchor projectile surface and the shear residual stress component acting on the projectile 
tangential to this surface are defined as confinement stresses. 

In the computation of contact/impact phenomena, more accurate contact stress data can 
sometimes be obtained from nodal force data which has been converted to contact stress by a 
tributary element method (Hughes, et al., 1976). Element stresses are averaged quantities at the 
center of the element and therefore are not strictly associated with the interface surface. 
However, nodal force data cannot be conveniently post-processed with TAURUS. Thus, only 
element stress data were considered in this analysis of confinement stress. Further, computed 
pressure contour data, although easy to interpret and invariant, was avoided in the analysis of 
confinement stress because pressure is an average of the three normal stress components. The 
difference in magnitude between confinement stress and pressure can therefore be substantial 
because of the influence of tensile stress which often occurs in the hoop direction. 

The following confinement stress components are of interest and pertain to any arbitrary 
plane containing the symmetry axis of the penetration model: normal stress (a x ), shear stress 
(r n ), both referred to the positive xz plane, for example, of the target as shown in Figure 26, 
as well as the maximum and minimum principal stresses (a^n and cr,^). The vertical normal 
stress components, o z and a y , and the maximum shear stress and the intermediate principal stress 
components are otherwise not required for the present discussion. 

The intersection of the vertical xz plane and the cavity is curved in the region of the 
projectile nose. The residual normal stress (a x ) and residual shear stress (r H ) components in this 
region are not strictly confinement stresses because they are not normal or tangential to the 
projectile. Local stress component data are actually needed in this region instead of the global 
stress component data being discussed. Local stress components are currently not conventionally 
processed by DYNA3D. Therefore, in the region of the projectile nose, the minimum principal 
stress (largest compressive stress) data should also be considered in the analysis of 
confinement stress in the absence of local stress data. 

An understanding of cavity expansion theory also helps to interpret numerically computed 
confinement stress data. From the results of this theory, similarity should be anticipated between 
the <r x and stress components. Cavity expansion theory assumes the strain energy in the 
target to be distributed exclusively in horizontal planes. In effect, this theory assumes a x to be 
identical to (in the vertical xz plane). Further, values for these components should tend to 
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Figure 26 

Stress block at center of a typical finite element in target model. 

be compressive. The extent to which calculated three-dimensional residual stress data differ from 
these observations, is the extent to which the essentially one-dimensional character of cavity 
expansion theory is in error. Conversely, cavity expansion theory suggests that stress 
components a x and should be dissimilar, and that they should be of opposite sign. This is 
because the theory effectively assumes that g ^,, is identical to the tangential or hoop stress 
component, which tends to be tensile. 

Discussion of Results 

Residual stress fields for the four components (g x , r^, g,^, g min ) are presented in 
contour data form and discussed for two anchor projectile simulations involving sandstone 
targets: the smaller fastener anchor projectile simulation and the larger anchor projectile 
simulation. For the respective input data files, see Appendix F, try ft 7 and run #1. Twenty 
DYNA3D simulations of penetration models were conducted during the course of this project. 
These two simulations were considered to be more realistic. Others were aimed at understanding 
the effects of various simulation parameters as previously described in the parameter study. 

Smaller Fastener Anchor Projectile. The residual stress results for the fastener anchor 
projectile problem are presented in Figure 27. The residual stress contour values are in units of 
kg/mm*s*(to convert to MPa, multiply the data in Figure 27 by 10‘ 3 ). Negative values in the 
figure refer to compressive stress. It is noted that, as required by rotational symmetry, the 
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Figure 27 

Residual stresses for smaller fastener anchor projectile in sandstone. 
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computed principal stress data shown are rotationally symmetric, while the normal and shear 
stress data shown are restricted to the xz plane. 

Irrespective of the particular component, the largest residual stresses occur adjacent to the 
embedded anchor projectile. Thus, confinement stresses are observed to be the largest residual 
stresses in the surrounding rock medium. Residual stress decays rapidly with distance from the 
embedded anchor projectile. Beyond two or three projectile diameters, the residual stress in the 
rock is negligible. This behavior is generally consistent with cavity expansion theory which also 
predicts the rate of decay of the residual stress. The particular cavity expansion theory presented 
in Appendix A predicts that residual stress decreases as 1/r 2 in the cracked region of brittle 
materials and as 1/r 3 in the elastic region beyond, where r is radial distance from the symmetry 
axis (see Figure A-4). 

In Figure 27(a), the maximum value for a x , the normal component of confinement stress, 
occurs halfway down the length of the embedded projectile shank (straight part). In Figure 
27(b), the maximum value for the tangential component of confinement stress, occurs 
slightly farther up on the embedded projectile shank. The contour pattern for the minimum 
principal stress component presented in Figure 27(c) is similar to the contour pattern for the 
normal stress component. The confinement stress diminishes almost an order magnitude from 
the maximum location, midway along the shank, to the anchor projectile tip. 

The maximum principal stress, shown in Figure 27(d), tends to be tensile, and 
approximates the hoop stress component (the a y component which is not shown) in the rock 
surrounding the projectile. As such, it is not a confinement stress per se, but it is indicative of 
the cracking potential of the rock surrounding the projectile, which can, in turn, affect 
confinement stress. The maximum value of this stress component occurs at the surface in the 
crater uplift area which is shown as a raised surface adjacent to the projectile. It diminishes to 
a depth midway along the cavity, whereupon it increases again toward the projectile nose. 
However, since cracking was excluded from these calculations, stress levels above the tensile 
strength of sandstone (0.1 to 4.5 MPa) are clearly unrealistic. 

The contour data presented above were calculated using a tunnel-to-projectile diameter 
ratio of 0.35. They are between the results for a ratio of 0.25, and the results for a ratio of 0.50 
(respectively, tries #9 and #8, Appendix F). The results for these two cases are shown in 
Figure 28, and were compared with the results in Figure 27(a). The size of the tunnel has little 
effect on the residual stress response if the ratio is near 0.25, the value arrived at on the basis 
of target penetration depth in the parameter study. 

The coefficient of friction influences the depth of anchor projectile penetration, but 
varying the coefficient from zero to 0.1 was found to have little influence on confinement stress. 
The residual stress field resulting from specifying no friction (try #12, Appendix F) is shown in 
Figure 29. It is very much the same as the residual stress field for a coefficient of friction of 
0.1 (Figure 27a). The difference is only in the location of maximum confinement stress which 
changes to a higher position on the shank of the anchor projectile due to greater depth of 
penetration. The anchor holding force should correspondingly increase due to greater 
penetration, but confinement stress (which is basically unit holding force) is basically independent 
of the coefficient of friction. 

In the parameter study, it was found that penetration depth was influenced by the 
hardening modulus of the rock when modeled with DYNA3D material model #3. In the present 
case, the hardening modulus was reduced from one percent of Young’s modulus to zero. The 
residual stress results are shown in Figure 30. When these data are compared with those for a 
hardening value of one percent (Figure 27a), a corresponding decrease in confinement stress of 
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approximately 25 percent is evident for the case of no hardening. The hardening modulus of 
rock, therefore, plays a significant role in the magnitude of confinement stress. 

Larger Anchor Projectile. The initial kinetic energy for the large anchor projectile 
simulation was two times greater than that for the fastener anchor projectile. However, the 
residual stress field, shown in Figure 31, when converted to a common measure,* was 
approximately the same as the residual stress field for the fastener problem (Figure 27). The 
similarity of these results is explained by the specification of the same material model. Both 
anchor penetration simulations employed the same elastic-plastic material model with the same 
yield strength and hardening modulus. That the initial kinetic energy of impact and projectile 
radius differed substantially, made little difference to the residual stress. Thus, residual stress 
essentially depends on material properties of the rock target. This is also observed in static 
cavity expansion theory, where the stress necessary to quasi-statically open a cavity of prescribed 
radius is the primary objective (see Figure A-2 in Appendix A). According to this theory, the 
stress is independent of dynamics, and primarily dependent on the strength of the target material. 
Penetration response of the projectile is determined separately by considering an auxiliary 
dynamic analysis of a penetration model. Thus, the use of appropriate material models for 
characterizing rock material response in simulation of penetration problems is very important 
according to the results from both numerical models and analytical models. 

The contour patterns for normal stress, a x , and minimum principal stress, a mm , are 
similar, as shown in Figures 31(a) and 31(c), respectively. This again indicates that the rock 
strain energy is primarily stratified and nearly contained in horizontal planes, as is assumed in 
cavity expansion theory. Alternatively, consideration of the vertical stress component, a z , may 
be ignored in favor of ct x (or a y ) as has been done in the present analysis of residual and 
confinement stress. The maximum value of the confinement stress component, a x , is located on 
the projectile nose. Slightly higher stresses occur just beneath the nose. The normal 
confinement stress component diminishes almost one order of magnitude along the projectile 
surface from the nose to the surface of the target. 

In Figure 31(b), the tangential or shear residual stress component, r n , is shown to be less 
in magnitude than the normal stress component by a factor of about six. The tangential 
confinement stress varies slowly along the projectile surface from a maximum value at the tip 
to a minimum value at the surface of the target. The negative contour values of r n at the 
interface indicate that the confinement stress acts upward on the target cavity and downward on 
the embedded projectile. 

Contour data for the maximum principal stress component, o max , shown in Figure 30(d) 
represent primarily a hoop stress component. This stress component is compressive (negative) 
along the lower half of the embedded projectile and tensile along the upper half. This indicates 
that radial cracking in the target would be limited to the upper portion of the projectile and 
certainly to the crater area at the surface. 


*To convert the labeled contour values in Figure 31 from g/cm*s 2 units to MPa units, multiply 
by 10’ 7 . 
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Figure 28 

The effect of tunnel-to-projectile diameter ratio on residual stress. 






Figure 29 

Normal stress, <r x , for zero coefficient of friction. 
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Figure 31 

Residual stress for larger anchor projectile in sandstone. 
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On a typical unit surface area of the projectile, the downward shear confinement stresses, 
Figure 31(b), are equilibrated by the upward component of normal confinement stress acting on 
the ogival nose. Figure 31(a). The free body diagram of the projectile, sketched in Figure 32, 
roughly illustrates the equilibrium of these two components of confinement stress. 



Figure 32 

Confinement stress distribution for larger anchor projectile embedded in sandstone. 


SUMMARY AND CONCLUSIONS 

Three-dimensional finite element modeling of penetration of ballistically embedded 
projectiles in rock was considered to establish a numerical procedure for investigating the 
influence of various parameters on the penetration process and on the resulting residual stress 
fields surrounding the embedded projectile. The explicit finite element program DYNA3D was 
employed and evaluated for these computations. Insight was also gained from a brief study of 
cavity expansion theories applied to projectile penetration of rock targets. In this preliminary 
study, the finite element penetration model was limited to normal projectile impacts and 
axisymmetric geometry. A deformable projectile model, nonlinear material behavior for the 
projectile and rock target, and a sliding interface with friction between the projectile and rock 
target were included in the penetration model. 
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The following conclusions pertain to construction of a penetration model for determining 
confinement stress around projectiles embedded in rock targets using DYNA3D. 

In general, an artificial tunnel must be included in the finite element model of the rock 
target to facilitate model construction and to minimize computer run time. The penetration 
depth, velocity, deceleration, and stress response were not affected by the tunnel if the tunnel-to- 
projectile diameter ratio was less than 0.25. In a similar manner, a small artificial bore at the 
center of the finite element model of the projectile may be used to enhance the numerical 
efficiency of the penetration model. 

The target may be divided into two parts to facilitate finite element subdivision. The 
inner part of the model must be finely subdivided to capture the high deformation gradients that 
occur near the target/projectile interface. The outer part of the target model may be coarsly 
subdivided, and its outer diameter may be limited to three to four times the projectile diameter 
without affecting accuracy, since stress fields decay rapidly with distance from the projectile. 

Specification of a sliding surface is essential in modeling the behavior of the 
projectile/target interface. Further, a tied sliding interface option was convenient for imposing 
displacement compatibility at the interface between the inner and outer parts of the target model. 

Twenty simulations of penetration models were conducted. The models included a smaller 
anchor projectile, representing a fastener, penetrating a sandstone target, and a larger rock anchor 
projectile also penetrating a sandstone target. All simulations were limited to normal projectile 
impacts. They were also limited to three-dimensional quarter models of axisymmetric penetration 
configurations to minimize simulation turnaround time using a Sun 4/260 workstation. 
Simulation times were found to easily allow for overnight turnaround time, and improved 
performance may be anticipated with newer workstation hardware. 

In summary the DYNA3D system proved remarkably facile in simulating complex 
penetration problems and calculating the structural dynamics of normally impacting, ballistically 
embedded projectiles and the associated confinement stress fields surrounding the embedded 
projectile. 

Direct attempts to validate the numerical penetration model using existing analytical 
solutions for penetration problems based on cavity expansion theory were not successful. The 
analytical solutions are generally too limited to serve as a benchmark. They cannot provide 
stress data results with sufficient detail commensurate with the detail in residual stress gradients 
surrounding the body of embedded projectiles. However, the results from a particular spherical 
cavity expansion theory were useful as a guide to interpretation of numerical results, and 
substantial qualitative agreement was found between numerical and analytical solutions for 
residual stresses surrounding embedded projectiles. 

In a further attempt at validation, an existing DYNA2D numerical solution of an 
axisymmetric penetration problem was replicated using the DYNA3D penetration model 
developed in this study. The DYNA3D prediction of penetration depth did approach the 
DYNA2D solution from below, as expected from the mathematical theory of the displacement- 
based finite element method. 

Parameter studies of the penetration models were conducted to assess the effects of 
alternative DYNA3D material models for rock targets including various target strengths, and to 
assess the effects of alternative coefficients of friction for a Coulomb model of the sliding 
interface between the projectile and rock medium. 

Of four DYNA3D material models indentified as potentially suitable for modeling rock 
targets, two were investigated in this study. The more successful model was an elastic-plastic 
strain hardening model. This model, however, does not allow for modeling the brittle behavior 
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of rock in tension. A second model, the soil and crushable foam model, does allow for a tension 
cutoff and brittle material modeling. However, this model failed to work when the tension cutoff 
value specified was other than zero. When it did work, a null solution was obtained for the 
residual stress state in the target material because the target cavity produced by the impact 
separated from the embedded projectile. 

The strain hardening modulus for the elastic-plastic material model for rock targets had 
the most influence on computed residual stress states. A value equal to 1 percent of Young’s 
modulus for the rock material, resulted in a 25 percent increase in confinement stress acting on 
the embedded projectile, and a similar decrease in penetration depth, when compared to the 
results assuming no strain hardening in an elastic-perfectly plastic material model characterization 
of the rock target. 

In contrast with the strength of the rock target material, the coefficient of friction at the 
projectile/target interface had only a minor affect on confinement stress. However, an increase 
in projectile deceleration and a decrease in penetration depth with increasing coefficient of 
friction did occur as expected. 


RECOMMENDATIONS 

Further studies to evaluate DYNA3D material models for rock targets are required. 
Material models must eventually account for cracking in rock media. Recently published 
DYNA3D documentation (Whirley and Hallquist, 1991) includes much improved definitions of 
material models and input parameters that will facilitate additional study of material models for 
rock. Material models accounting for more realistic behavior including anisotropic in-situ 
fracture conditions which could greatly influence the penetration process may yet have to be 
developed. The experience of the present study is that rock target material modeling remains 
the most important issue in successful numerical simulation of rock anchor penetration. 

Further attempts at validating the numerical penetration model should be pursued and 
should incorporate the experimental data compiled by Goldsmith and colleagues under NCEL 
sponsorship, Brunette and Goldsmith (1990), Goldsmith and Zhou (1989), and Zhou and 
Goldsmith (1990). 

Additional studies are needed with larger three-dimensional models that do not exploit two 
planes of symmetry for convenience. Supercomputers should be used. DYNA3D exists on 
many of these platforms, and penetration models developed initially on workstations can be easily 
moved to supercomputing environments where performance increases of two orders of magnitude 
may be anticipated. 

Parameter studies are needed to evaluate the effect of finer mesh subdivisions on the 
dynamic response of projectiles and on confinement stress. Finer model subdivisions nearer the 
projectile would improve the definition and accuracy of confinement stress predictions. 

There are some post-processing limitations with TAURUS. Local stresses as well as 
nodal point force data should be made available for assessing confinement stresses more 
accurately. 

Upon refinement of the finite element penetration model as suggested in the 
aforementioned recommendations, further parameter studies should be resumed to investigate the 
relationship between residual stress fields and the properties of seafloor rock and PEA projectiles. 
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Appendix A 

A SPHERICAL CAVITY EXPANSION THEORY 

An overview of a spherical cavity expansion theory applied to projectile penetration is 
presented. The theory facilitates understanding of fundamental behavior in rock anchor projectile 
problems, and is used in the present project as a guide to interpreting results from the numerical 
models. 

Introduction 


The basic idea is to analyze the penetration of an object in a semi-infinite medium by 
simulating it as a cavity expanding in an infinite medium. 

As an example, consider that the medium is elastic-perfectly plastic, incompressible, and 
obeys the von-Mises yield condition, and that quasi-static expansion of a cavity from radius 1 = 
0 to r = a is to take place under the action of an unknown stress o s as shown in Figure A-l. 
The stress is to be determined by cavity expansion theory. 

In this case, cavity expansion theory yields the following radial stress at the cavity radius, 
r = a (see Schulz, 1978, or Bernard and Creighton, 1976), 

o, = |°y 1 + ln(| e y) (Ala) 


where a y and « y are the unaxial stress and strain, respectively, at the yield point. The maximum 
shear stress in this material is assumed to be 



(A-lb) 


An auxiliary analysis of a simple penetration model is developed next. Assume during 
penetration that the stresses in Equation A-l act over the hemispherical tip of a rigid spherical¬ 
nosed projectile. The total drag force, F s , acting on the projectile can then be obtained by 
integration over the hemisphere, as follows. 



Conservation of momentum for the projectile requires that 

F f = m v (A-3) 


where m is the projectile mass and v is the projectile velocity, and a superscript dot indicates 
differentiation with respect to time. 


A-l 





Figure A-l 

One-dimensional coordinate system, stress and material behavior. 


Equations A-l, A-2, and A-3 can be combined to yield a first order differential equation 
for the velocity as a function of time, which can easily be integrated. A second integration will 
yield the penetration depth, z, of the projectile as a function of time. The results are: 


v(t) = 


na 2 ( tc \ 


t - v n 


(A-4) 


z(t) = 


n a 2 ( it \ t 2 


(A-5) 


where v 0 is the initial impact velocity at t = 0. 

The crucial part of the penetration model analysis is knowing the radial stress on the 
cavity surface. Once known, the depth of penetration can then be determined similarly as shown 
by the above simple penetration model example. 

Spherical Cavity Expansion Theory for Infinite Brittle Medium 

A spherical cavity theory is presented for determining the radial stress at the cavity radius 
for brittle materials such as concrete and rock. The presentation is based on a re-derivation of 
the equations published by Forrestal and Longcope (1990) with extension of the results to plain 


A-2 





concrete. The equations presented in this brief publication were confirmed, and their re¬ 
derivation in this study served to further explicate the theory. 

Referring to Figure A-2, the problem addressed is to find an analytical expression for the 
radial stress a r (a). The following primary assumptions are made: 

• Spherical symmetry (<r* = a e ) 

• Homogeneous and isotropic material 

• Quasi-static cavity expansion 



Material descriptions for the three regions indicated are required. The elastic region is 
characterized by the following constants: 

K = bulk modulus 
E = Young’s modulus 
¥ — Poisson’s ratio 
T = tensile strength 

The cracked region is characterized as follows. Radial cracks are assumed to exist and, 
therefore, = <r t = 0. Otherwise, the region is elastic with the same parameters as the elastic 
region. It is this region that specializes the theory for brittled media. 


A-3 






The plastic region is characterized by the following constants. 

Y = compressive strength 

X = pressure-dependent constant for shear strength (0 £ X «s 1) 
r = shear strength which depends on confinement or mean stress, and is calculated by 

t = (3-X)Y/3 (A-6) 


Elastic Region, r s R. The boundary value problem describing the elastic region is: 

Governing equation of equilibrium: —- + — (o r - o 0 ) = 0 (A-7a) 

dr r 


Boundary condition: o t (R) s p 


(A-7b) 


where the radial stress P will be related to the tensile strength of the elastic material on the verge 
of cracking. 

The solution of the boundary value problem, which is given by Hill (1950), is: 


Radial stress: o r 



(A-8a) 


Circumferential (tension) stress: o. 


PR 3 

2r 3 


Meridional stress: = o 0 


Radial displacement: u = 


1 + v PR 3 


2 E 


At the boundary r = R, the stresses are evaluated as 


°r - P 


and 


= = 



(A-8b) 

(A-8c) 


(A-8d) 


(A-9a) 


(A-9b) 
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(A-10) 


Further, due to cracking of the elastic material at the boundary, 

- - T 

Therefore, the radial stress P at the boundary can be determined. 

p = 2T (A-11) 


In summary, for the elastic region, the radial stress at r = R is twice the tensile strength, 
T, and diminishes as 1/r 3 as shown in Figure A-3. 



Figure A-3 

Radial stress in elastic region. 


Cracked Region, c s r s: R. The boundary value problem describing the cracked region 


do 2 

Governing equation of equilibrium: — + — o r = 0 

dr r 


(A-12a) 


Boundary condition: o r (R) = 2T 


(A-12b) 
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The solution to the boundary value problem, by the method of separation of variables, 


is: 


a r = 


2 TR 2 


(A-13) 


Apply the boundary condition at r = c, that the radial stress equals the compressive 
strength, Y. 


o t (c) = Y 


(A-14) 


By eliminating T in terms of Y in Equation A-13, an alternative expression for radial stress is 
found. 



(A-15) 


Equating the two expressions in Equations A-13 and A-15 for radial stress, the following 
relationship is obtained. 



(A-16) 


The following remarks stem from Equation A-16. 

• The ratio of radii defining the limits of the cracked region is fixed by the 
tensile and compressive strengths of rock. 

• It may be supposed that the definition of a brittle material becomes one in 
which the compressive strength is more than twice the tensile strength so that 
the radii ratio does not exceed one. 


• The values of c and R are obtained after consideration of the plastic region and 
completing the solution of the problem. 

To obtain the radial displacement for the cracked region, the pressure-volume relation is 

applied 

p = Ke (A-17) 


Here, p is the mean stress which, for the cracked region, is, 

1 

p = — o, 

3 r 


(A-18) 
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(A-19) 


and «is the percent volume change which, for the cracked region, is, 

e = e t + 2e e 


The strain-displacement relations are, 


e. = - 


du 

dr 


(A-20a) 


and, 


e„ = - - 


(A-20b) 


Therefore, using Equations A-17 through A-20, the following relation is obtained: 

h = + (A-21) 

3 [dr r] 

When Equations A-15 and A-21 are combined, the following equation governing the radial 
displacement is obtained: 


Pressure-volume relation in terms of u: 


du 2u 
dr + r 


Yc 2 

3Kr 2 


(A-22a) 


Boundary condition from elastic region: u(R) 


•(W) 


PR 


(A-22b) 


The solution of this boundary value problem, by the method of integrating factors, is: 


u = 


(1 » v)TR 2 + Yc 3 


3K 


Yc 2 1 


3K r 


— (csrsR) 


(A-23) 


This completes the analysis of the cracked region. The stress distribution for the cracked 
and elastic regions is shown in Figure A-4. In the cracked region, the radial stress varies as 
1/r 2 , compared to 1/r 3 for the elastic region. This shows that the presence of radial cracks 
impedes the decay of radial stress with distance from the source in an elastic medium. 

Plastic Region, a £ r s c. The governing equation of equilibrium for the plastic region 
is 


do, 2(o r - o fl ) 
dr r 


(A-24) 


A-7 






= radial stress 
= radial axis 
= point at cracked region 
= tensile strength 
= point at elastic region 
= compressive strength 



Figure A-4 

Radial stress in cracked and elastic regions. 


The material is in a state of plastic deformation. Accordingly, stresses must obey the pressure- 
dependent failure surface relation shown in Figure A-5. This yield surface is described by 


(o t -o e ) = A-p + t 


(A-25a) 


where 


P “ -(o r + 2o e) 


(A-25b) 


Eliminating a 6 in terms of a t , the governing equilibrium equation becomes, 

do of (3 o r - 2 t \ 

+ - x \—--Ut = o 

dr r[ ( 3+2A ) 

To simplify the algebra, a new constant a is defined, 


a = 


6 

3 + 2 A 


(A-26) 
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Figure A-5 

Constants model for pressure-dependent shear strength. 


Thus, the boundary value problem for the plastic region is: 


Governing equation of equilibrium: 



tt T 

r 


(A-27a) 


Boundary condition: °r<c) = Y 


(A-27b) 


The solution of the boundary value problem, by the method of integrating factors, is: 


Y 

* T 


i r ■ (¥)] 


(A-28) 


This expression for radial stress in the plastic region when evaluated at r = a gives the 
required solution: 


o s h o r (a) 


(A-29) 


However, the ratio c/a remains unknown. The required solution for the ratio c/a is outlined in 
the following steps: 
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(1) Invoke the mass conservation equation in Eulerian coordinates and the displacement 
boundary condition at r = c from the cracked region. 


d 3 ot 2 

Mass conservation: — [(r - u) 3 ] = —-— 

dr Po 


(A-30a) 


Boundary condition: u(c) = — + (R - c) 

7 Ec 2 3K 


(A-30b) 


(2) Approximate the mass density ratio for small volume changes as follows, and 
express e in terms of the radial stress, 


P i i i 

— = - « 1 + e = 1 + 

Po 1 " e 


11 

i 

cn 

a 

O ] 

!m 

rj 3 X j 


(A-31) 


(3) Substitute Equation A-31 into Equation A-30a and integrate both sides from a to 
c: 


[c -u(c)] 3 - (c 3 - a 3 ) 


{ 3k ) K 


3 Y 


X K(3 - X a) 


-ID" 


(A-32) 


(4) Expand the left-hand side of Equation A-32, neglecting second and higher order 
terms, 

[c-u(c)] 3 « c 3 - 3 c 2 [u(c)] (A-33) 

(5) The resulting expressions for the required ratio, after some algebra, are, 

First order solution: j * 6 = — (1 - v) f —— 1 — (A-34a) 


\3 a ( Y \ lf2 Y 

* p I 


Second order solution 


- (!) * 


3 n ■ 

k E 


(3 - X) - (3 + 2X) 


l V* 2 ' 3 


(1) 

IpJ 


(A-34b) 


In these expressions, T is the tensile strength, E is Young’s modulus, and v is Poisson’s ratio. 
This completes the solution of the plastic region and the spherical cavity problem. 

Summary. The pressure necessary to open a spherical cavity from zero initial radius to 
radius r = a is, 


o, s o t ( a) = 


ter-ml 


(A-35) 
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the first order solution 


the uniaxial compressive strength 

the pressure dependent shear strength parameter (0 £ 1 £ 1) 



Appendix B 

ANALYSIS OF A DYNAMIC CYLINDRICAL CAVITY EXPANSION THEORY 


The penetration theory presented by Zhou (1988) was analyzed in an attempt to calculate 
depths of penetration. The units used in this analysis are kilogram, meter, and seconds. The 
equation numbers are Zhou’s. The variables are defined as follows: 


a : radius of projectile 

R : radius of nose 

A : area of cylindrical section 

Vo : initial velocity 

Cl : quasi-static coefficient 

Y : yield strength 

C2 : inertial coefficient 

Y1 : effective yield strength 

d : diameter of projectile 

Zf : depth of penetration 

E : Young’s modulus 

a2 : guess 

Et : strain hardening 

0 : material constant 

g : gravity 

«1 : strain in locked elastic region 

1 : projectile length 

e2 : strain in locked plastic region 

L : length of nosecone radius section 

0 : ogival half-nose angle 

m : mass of projectile 

pi : density of elastic region 

n : number of terms in summation 

p2 : density of plastic region 

percent po : assume a % porosity for pi 

po : initial density of target 

Ps : shear resistance 

ps : density of steel projectile 


The initial values given the variables are as follows: 
d := .004 R := .010 1 := .035 po := 2130 Vo := 100 ps := 7820 

«2 := .01 percentpo := .18 n := 1 ..10 E : = 13.6210 9 Y := 4510 6 


by the following relationships: 
pl 

1 - a2 
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The densities of the elastic and plastic regions are given 
pi := po (l +percent_po) p2 : = 








The relationship used for the effective yield strength, strain hardening, and the strain in 
the locked elastic region are as follows, respectively: 


Y1 := Y • 1 - — Et : = 1.1 E el := — 
E E 


(A. 20) 


The material constant, 0, as defined by Zhou: 

* . 3•Y [el 


The ogival half-nose angle, 6, is developed: 

a : = — L : = ^R 2 - (R - a) 2 0 : = — - a sin ——— 

2 V 2 R 


The quasi-static term and the initial coefficients ds defined by Zhou: 

Cl := — L a-cos(0) - L + —— - + -sin(0)-cos(0)l 
cos(0) sin(0) 2 2 


(2.42) 


2L f _ L 
a 2 (cos(0)) 2 sin(0)_ 


po L[ sin(6) - cos(6J » 1 1.^ 
2 sin(0) + cos(0) - 1 


(2.43) 


_ ln(a2) • tan(0) - * - sin(0) 

2 L • cos(0) 


-- — — - — sin(0) • cos(0) ln(a2) 

a 2 • sin(0) • cos(0) 2 2 2 2 2 


■ |l - -isinO)] • l n[ s ‘ p(e )- cos(e ) tl l-[- + Ism(0)cos(6) 
L J sin(0) + cos(0) -1 2 2 


C2 := C2a + C2b (Sign may be + or -; it was not given.) 
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The mass of the projectile, m, is given as follows: 
A:=it-a 2 m:=[(l-L)-A + 


(2.20a) 

(2.41) 

The values obtained for the variables are as follows: 
m = 0.003 0(180/*) = 36.87 

po = 2130 pi = 2513.4 p2 = 2538.788 «2 = 0.01 

Cl = 0 C2a = -3720.066 C2b = 2099.201 

Ps = 5.829 • 10 17 Zf = 4.92145418 • 10' 13 

This is an incorrect Zf depth. The shear resistance Ps appears to be unusually large. 
This leads to a check of the shear resistance term or more specifically, the summation of the 
shear resistance term. 

The shear resistance summation term was computed for several values of a2. This term 
is computed as follows: 

«2 := .1, .2 .. .5 n := 1 .. 25 





-LA -ps 
4 
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a =. 1 

a = .2 

a =.3 

a = .4 

a =.5 

0 

0 

0 

0 

0 

- 1.081 

- 0.618 

- 0.413 

- 0.291 

- 0.207 

2.279 

0.921 

0.503 

0.306 

0.195 

- 4.821 

- 1.354 

- 0.603 

- 0.321 

- 0.186 

10.718 

2.081 

0.756 

0.353 

0.186 

- 24.911 

- 3.345 

- 0.992 

- 0.405 

- 0.194 

60.006 

5.594 

1.354 

0.483 

0.21 

- 148.768 

- 9.665 

- 1.91 

- 0.595 

- 0.234 

377.636 

17.158 

2.768 

0.75 

0.267 

- 977.648 

- 31.157 

- 4.104 

- 0.967 

- 0.31 

2573.52 

57.656 

6.2 

1.268 

0.366 

- 6871.898 

108.412 

- 9.518 

- 1.688 

- 0.437 

18578.048 

206.638 

14.813 

2.278 

0.53 

- 50771.228 

- 398.493 

- 23.323 

- 3.109 

- 0.648 

140076.327 

776.301 

37.098 

4.285 

0.8 

- 389733.04 

- 1525.759 

- 59.533 

- 5.957 

- 0.996 

1.093 10 6 

3022.237 

96.284 

8.345 

1.249 

- 3.083 • 10 6 

- 6028.02 

- 156.803 

- 11.771 

- 1.577 

8.753 • 10 6 

12097.682 

256.943 

16.705 

2.003 

- 2.499 • 10 7 

- 24413.887 

- 423.375 

- 23.839 

- 2.557 

7 . 168 - 10 7 

49515.906 

701.111 

34.19 

3.282 

- 2 . 066 - 10 ® 

- 100884.335 

- 1166.326 

- 49.258 

- 4.229 

5.977 • 10 ® 

206394.848 

1948.271 

71.26 

5.473 

- 1.736 10 9 

- 423854.969 

- 3266.798 

- 103.479 

- 7.109 

5.059 10 9 

873463.732 

5496.727 

150.788 

9.267 


This shows the shear resistance term diverges. This numerical result prompted a theoretical 
check of the series term for shear resistance using the n th term test for a divergent series. 


The n th -Term Test for a Divergent Series. The following series is a term in the 
formula for target resistance in the cylindrical cavity expansion theory presented by Zhou (1988). 



(B-la) 


where, 

a = l_fi (0 < a < 1) (B-lb) 

P 2 . 

The numerical evidence indicates this term to be divergent. Analyzing it by the n th -term 
test for divergence, it is first simplified by choosing a specific value of a (with no loss of 
generality). A convenient value is, 
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a 


(B-2) 


4 


The series simplifies to, 


E^-rC-ir 1 < B - 3 > 

n-1 n 

Allowing the n th term to pass to infinity, 

lim (_!)“-! = lim — (-1)®' 1 (B-4) 

n-°® n 2 n - «> n 2 


Applying L’HospitaTs rule, 


n — n 2 


lim 2 nl n(-l) n - 1 + 2°(n -1)(- l) n ~ 2 
n-« 2 n 


(B-5) 


Applying it again, 


lim —— (— l) n ' * 
■— n 2 


lim 2 n - 2 n(n-l)(-l) n l + 2n°' 1 (n-l)(-l) p ~ 2 + ... 

n-« 2 


(B-6) 


This limit oscillates between + and - infinity, and therefore, 

limZ-ll * 0 (B-7) 

n -- n 2 


Therefore, the series B-3 diverges. By inspection, this result is typical of any value of a in its 
range (0 to 1). Therefore, the series B-l is divergent by the n th -term test. 
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Appendix C 


OVERVIEW COMMENTS ON DYNA3D 

Three features of the computer program DYNA3D are discussed which are relevant to 
understanding its application to penetration problems involving anchor projectiles in rock. These 
features relate to the program’s capability to handle highly transient events involving material 
and geometric nonlinearity and interface sliding capability in a robust manner. The discussion 
equally applies to the two-dimensional program, DYNA2D. Further, the discussion given is 
introductory, and is generally meant for potential users of the program who are only somewhat 
familiar with finite element theory and jargon. A detailed description of DYNA3D is provided 
in the user manual (Whirley and Hallquist, 1991). 

Explicit Time Integration 

DYNA3D is an explicit finite element computer program for three-dimensional nonlinear 
dynamic problems. The word explicit refers to the type of direct step-by-step method of 
numerically integrating or discretizing the system of nonlinear equations of motion. These 
equations are second order differential equations in time. Specifically, DYNA3D employs the 
explicit central difference algorithm in its solution. The central difference algorithm advances 
to the next time step using only known information calculated in the previous two time steps, and 
because it relies solely on known information while advancing the solution, the algorithm is 
termed explicit. 

In contrast, an implicit finite element computer program requires not only information 
calculated previously, but also requires unknown information at the end of the time step to 
advance the solution. To do this it must also solve a system of simultaneous algebraic equations 
which result from the implicit algorithm’s discretization of the equations of motion. This system 
is often large, and it must be solved each time step because of the nonlinearity present in the 
equations of motion. 

There are, however, two advantages to implicit algorithms. The time steps may be 
substantially larger than for explicit algorithms, and the numerical process may be 
unconditionally stable, whereas the explicit algorithms are only conditionally stable. The latter 
advantage strictly applies to linear problems, however. 

A characteristic of explicit temporal integration algorithms is that they must employ very 
small time steps to maintain stability of the numerical process. This is a consequence of their 
conditional stability. The time step size is based on a characteristic length of the smallest 
element in the finite element model. DYNA3D monitors the change in the characteristic length 
during the solution and automatically adjusts the time step size so that stability is maintained. 
The user has no control of how small the time step becomes during the solution, except of course 
initially in the construction of the finite element model. 

The net result of employing the explicit time integration algorithm is that DYNA3D can 
handle highly transient events and complex problems in a fairly straightforward, if not brute 
force, manner. Abrupt changes in nonlinear material behavior, deformation and loading, all of 
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which are reflected in the current state of the equations of motion, are generally captured 
satisfactorily by small time steps and explicit methods. However, very large numbers of time 
steps are required even for simulations of highly transient events. 

Constant Stress Finite Elements 

In addition to the straightforward simplicity of an explicit integration algorithm for 
solving complex nonlinear equations of motion, DYNA3D employs finite elements which allow 
for simple evaluation of integrals in the finite element formulation, such as the internal force 
vector integral. This is motivated by the need to minimize the compute time within each time 
step since large numbers of time steps are required. To this end, the DYNA programs employ 
constant stress finite elements for solving continuum problems; specifically, a four-node 
quadrilateral element in DYNA2D, and an eight-node solid element in DYNA3D are used. 
These are very basic displacement-based finite element formulations. 

Further, DYNA3D defaults to a one-point Gaussian quadrature rule at the center of each 
element to numerically integrate the required integrals. The required integrals are therefore 
under-integrated for the sake of minimizing the compute time within each time step. The 
material constitutive equations are therefore only evaluated at the center of each element. Also, 
a form of numerical instability called "hour glassing" becomes a real concern due to the fact that 
element stiffness matrices are not fully evaluated (they are rank deficient). They may be prone 
to assuming spurious shapes, actually resembling an hour glass in some cases. This instability 
is generally controlled in DYNA3D and is otherwise not a practical concern to the user. On the 
other hand, specification of full quadrature rules generally results in prohibitive overall compute 
times for problem solution and is often not a desirable option. 

Special Handling Techniques for Contact/lmpact 

Originally, the DYNA3D algorithm for integrating the equations of motion for nodes 
which either impact or release on interface surfaces was based on a contact/impact algorithm for 
implicit temporal integration developed by Hughes, et al. (1976). Special equations for velocities 
and acceleratio s in the basic temporal integration algorithm were developed to manage the jump 
conditions at impact and release of sliding interface nodes. In simple planar impact 
demonstration problems, they were shown to dramatically improve computational results by 
eliminating ringing in the impact stress data. Other demonstration problems were subsequently 
solved in two dimensions involving Hertz contact of curved surfaces, and large deformation 
contact problems of curved surfaces. These equations were originally incorporated in DYNA3D 
(Hallquist, et al., 1976) after first adapting them to the explicit temporal integration algorithm. 

It was reported (Hallquist, et al., 1977) that in hydrodynamic applications, where the 
material in constant stress elements support large confinement pressures well above the material 
strength, "hourglass" or "zero-energy" mode shapes contaminated the solution and was a 
recurring problem. It was not clear as to whether the aforementioned special handling techniques 
were the cause of the problem. However, it was eventually concluded in the evolution of the 
DYNA3D contact/impact capability, that for general application these special handling techniques 
were too expensive and they have since been dropped. 

To solve contact/impact problems, the following procedure was developed to handle 
sliding surface interfaces in DYNA3D. For each sliding interface prescribed in the finite element 
model, DYNA3D runs through the following non-iterative update steps for every time step: 
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1. Slave node search 

2. Slave element search 

3. Tensile interface force calculation 

4. Add slave mass to master mass 

5. Momentum calculation 

6. Effect of slave element stress on master surface acceleration 

7. Effect of slave mass on master surface acceleration 

8. Slave node accelerations and velocities 

9. Time step size necessary to prevent violation of impenetrability 

The ability to prescribe various sliding surfaces in DYNA3D makes it particularly 
applicable to the problem of penetration of a target by a deformable anchor projectile. In 
particular, the following types were employed in the present project: the "tied" (type 2) slide 
surface for joining two parts of the finite element model together to saftisfy displacement 
compatibility at the interface without concern for different mesh subdivision; and the "sliding 
with separation and friction" (type 3) slide surface for the interface between the anchor projectile 
and the rock medium. In the latter case, a penalty formulation is employed to manage the 
implementation of impact/release of adjacent surfaces, and a Coulomb friction model is employed 
for tangential or sliding motion of the adjacent surfaces. When a slave node penetrates an 
element on the master surface thereby violating the impenetrability condition, a "penalty" force 
is applied to the node so that it satisfies the constraint. 
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try #18, friction = . 1 
dn3d 

Id 1 lp 2 .2 0 .2 35. 

Id 2 Ip 2 0 .7 .6 .7 
Id 5 Ip 2 2. 35. 2. 5.864 
lap . 1 0 -8. 5.864 
Id 9 lp 2 .01 35 2. 35. 

Id 25 lp 2 .02 0 150. 0 
Id 27 lp 2 0 -150. 160 -150. 
Id 22 lp 2 .7 0 .7 -160. 

Id 24 lp 2 15. 0 15. -160. 

Id 28 lp 2 150 0 150 -160. 
Si 2 sv fric 0. ; si 3 tied; 


plane 2 0 0 0 0 1 0 .001 symm 

0 0 0 1 0 0 .001 symm 

term .5e-3 

plti 2.50e-5 

mat 1 3 head 

steel case - 6dac 

ro 7.8e-6 

e 209e6 

pr .3 

sigy 1.2e6 
etan l.e6 
mat 2 3 head 
sandstone 
ro 2.13e-6 
e 13.62e6 
pr .2 

sigy ,02e6 
etan 0 
endmat 

part 2 5 9 1 1 2 24 

drag rota 5 0 0 0 0 0 1 90 
si 2 1 1 2 2 2 2 m 
velo 0 0 -100e3. 
end 

part 22 27 24 25 2 102 4 
drag rota 5 0 0 0 0 0 1 90 
si 1 1 1 2 2 1 2 s 
si 1 1 2 2 2 2 3 s 
b 2 1 1222111111 
end 

part 28 25 24 27 2 34 -12 2.0 

drag rota 5 0 0 0 0 0 1 90 

si 1 1 2 2 2 2 3 m 

b 1 12 12 1 111111 

end 

end 


Problem identity statement. 
Specifies a DYNA3D input tile 
Starts the line definitions for the 
projectile geometry. Lines 1, 2, 
5, and 14 describe the projectile 
geometry. 

Lines 22, 24, 25, 27, and 28 
describe the two part target 
geometry. 


Defines slide interface #2 as a frictional sliding 
surface with voids (gaps) permitted. Slide interface 3 
ties the two parts of the target together. 

Planes 1 and 2 are defined to 
exploit the two planes of symmetry. 

Simulation stop time. 

Simulation plot interval. 

Defines material #1 as DYNA3D material model 3. 
Header comment as indicated in previous line. 

Density 

Young’s Modulus 
Poisson’s ratio 
Yield stress 
Hardening modulus 

Defines material #2 in the same manner as material 
# 1 . 


Ends material definitions. 

Defines a 2D projectile part of four lines, assigns 
material #1, and makes a 2 by 24 mesh. 

Creates a 5 piece solid by rotating the mesh 90 degrees 
Positions sliding interface #2 and makes it the master. 
Initial velocity for the part in z direction. 

Ends part definition. 

Creates the fine 2D target mesh. 

Creates a solid from the 2D mesh. 

Locates slide interface #2 and defines it as slave. 
Locates slide interface #3 and defines it as slave. 
Constrains the nodes at the bottom of the target. 

Ends part definition. 

Defines the coarse 2D target mesh. 

Creates a solid from the 2D mesh. 

Position the master for slide interface #3. 

Constrains the nodes at the bottom of the target. 

Ends part definition 
Ends input file 
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The following statements describe a sample session for creating a simulation run for 
DYNA3D. The sample problem is the large anchor projectile. The first action is to create and 
save the following lines in a ASCII text editor. The lines that begin with a c are comment lines 
and are not required for the program to run. 


Run #1, Large Anchor Projectile Baseline 

c friction = . I 

c target material model = 3 

c target hole radius = .3175 

c units are g, cm, sec. 

dn3d 

c states that a DYNA3D input file will be created 

c next 5 lines describe 1/4 of the projectile profile 

Id 1 Ip 2 .254 0 .254 17.78 

Id 2 lp 2 0 .508 1.27 .508 

Id 3 lp 2 1.27 17.78 1.27 4.211 

lap 0 0 -6.35 4.211 

Id 4 lp 2 0 17.78 1.27 17.78 

c next 5 lines describe 1/4 of the target profile 

Id 21 lp 2 0. 0. 25.4 0. 

Id 22 lp 2 .3175 0 .3175 -63.5 
Id 23 lp 2 0. -63.5 25.4 -63.5 
Id 24 lp 2 10.16 0 10.16 -63.5 
Id 26 'p 2 25.4 0. 25.4 -63.5 
c Define the sliding interfaces 
si 1 iv fric . 1; si 2 tied; 

c Sliding interface 1 is sliding with voids and has . 1 coef. of friction 

c Sliding interface 2 is a tied boundary 

plane 2 0 0 0 0 1 0 .001 symm 0 0 0 1 0 0 .001 symm 

c describes the planes of symmetry 

c The termination time and plot interval are set. 

term 5e-3 plti le-5 

c Start define material #1 

mat 1 3 head 

4340 steel projectile 

c assign material model #3 and puts in a heading statement 
c density = ro = g/cm A 3 
ro 7.8 

c Young’s modulus = e = g/(cm*s A 2) 
e 2.00el2 

c Poison’s ratio = pr 
pr .3 

c Yield stress = sigy = g/(cm*s A 2) 
sigy 17.9e9 

c Tangential modulus = etan = g/(cm*s A 2) 

etan 10.0e9 

mat 2 3 head 

sandstone 

ro 2.16 

e 44e9 

pr .3 

sigy ,234e9 
etan .440e9 
endmat 
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part 2 3 4 1 1 2 14 
drag rota 5 0 0 0 0 0 1 90; 
si 2 1 1 2 2 2 1 m 

c initial velocity of the projectile in cm/sec. 

velo 0 0-10000. 

end 

part 22 23 24 21 2 100 4 
drag rota 5 0 0 0 0 0 1 90; 
si 1 1 1 2 2 1 1 s 
si 1 1 2 2 2 2 2 s 
b 2 1 1222 111111 

c Constrain the boundary in the x, y, and z directions 
end 

part 26 21 24 23 2 25 -3 2.0 

c defines a four sided region, lines 26, 21, 24, 23, assigns material 2 
c 34 elements between lines 21 & 23, 12 elements between lines 24 & 26 
c -12 (- sign) & the 2 defines the ratio of the first to the last 
c segment between lines 24 and 26. 
drag rota 5 0 0 0 0 0 1 90; 
si 1 1 2 2 2 2 2 m 
b 1 12 12 1 111111 

c Constrain the nodes at the bottom boundary in the x, y, 

c and z directions. 

end 

end 

The above file was saved as a file named runl for this example. 

Next, the INGRID program is executed by taping ingrid at the prompt. (The program 
is case sensitive; note the type case when the file names are typed.) The INGRID program then 
asks for an input and output file name. In this example, respond with: 


i=runl o=runlo 


INGRID then loads the ASCII file runl and prompts the user for the graphics device 
driver. Select number 2 for screen output on the Sun IV in graphics window 2. At this point, 
any of the interactive commands can be entered. To create an input file for DYNA3D type: 

continue 

This causes INGRID to create a DYNA3D input file named the output file name given 
above. In this example, the file name is runlo. INGRID then exits. To start the DYNA3D 
program type: 

dyna3d 

DYNA3D responds with a prompt for input and output file names. In this example, the proper 
response is: 

i=runlo g=runlp 
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The input file name (i = ) is the same as the output file name from INGRID. The output file 
name (g=) is runlp and contains all the graphics information for the simulation. DYNA3D 
terminates when processing is finished. To view the results from the DYNA3D simulation, the 
post-processor TAURUS is used. To start TAURUS, enter the following command at the 
prompt in graphics window 2: 

taurus 

TAURUS prompts for the input file name. The proper response for this example is: 
g=runlp 

TAURUS then prompts for the graphics driver. Enter 2 for a Sunview driver (1152 x 900 
resolution). TAURUS is now ready for interactive commands. The following commands created 
some of the plots shown in this report for the large anchor projectile. 

state 500 

angle 1 ry 75 rx 30 xscale 5 yscale 5 zscale 5 
xtrans -15 ytrans 0 v 
m 2 

contour 8 
contour 9 
restore 
state 500 
rx 90 rz 180 

xscale 5 yscale 5 zscale 5 

angle 1 

ytrans -20 

xtrans -5 v 

m 2 v 

contour 8 

contour 9 

rx 90 

ytrans 5 

contour 9 

contour 8 

phs2 mat! 1 1 gather 
print mtime 3 1 1 
mtime 6 1 1 
mtime 9 1 1 
phsl 
exit 


This concludes the example problem. To create a file that can be printed at a postscript 
printer, create the lines above in a ASCII editor and add a 4 to the beginning of the file. This 
causes TAURUS to send all graphic images to a file in postscript printer format. Save the file 
under the filename lpplot. The response to the TAURUS prompt for input and output file names 
now becomes: 
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g=runlp c=lpplot 


This executes all the commands in the file lpplot and creates a file of the images that would 
otherwise be printed to the screen. TAURUS responds with the file name that the images are 
stored in. The file name may be dig.ps or dig.psa or dig.psb or etc. If the file is named dig.ps, 
the file can be sent to the printer by entering the following line at a command prompt once 
TAURUS is done processing: 

lpr dig.ps 

This causes the line printer to start printing. It will take a few minutes to print the file. The file 
used to create some of the plots for the Zhou penetration problem is given below: 

4 

state 500 

angle 1 ry 75 rx 30 xscale 5 yscale 5 zscale 5 
xtrans -35 ytrans -20 v 
m 2 

contour 8 
contour 9 
restore 
state 500 
rx 90 rz 180 

xscale 5 yscale 5 zscale 5 

angle 1 

ytrans -45 

xtrans -60 v 

m 2 v 

contour 8 

contour 9 

rx 90 

ytrans 60 

contour 9 

contour 8 

phs2 matl 1 1 gather 
print mtime 3 1 1 
mtime 6 1 1 
mtime 9 1 1 
phsl 
exit 
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The following is a list of the INGRID input files. For Zhou’s penetration problem, each 
file was called a try. The first four tries for the simulation concentrated on developing the 
slideline technique. Try #5 contained a large plot interval, and, as a result, the response of the 
projectile could not be tracked. The boundary in try #6 was not constrained. It was not until 
try #7 that the results were acceptable. Only the input data deviations from try #7 are 
highlighted. 

try #7: Friction = .1, Target hole radius = .7, etan = 1%, boundary constrained 
try #8: Target hole radius = 1 mm 
try #9: Target hole radius = .51 mm 

try # 10: Switch master and slave designation on slide boundary 2 

try #11: Friction on slide boundary 2= .05 

try tt 12: Friction on slide boundary 2=0 

try #13: Target material model #5. pc=0 

try #14: Target material model #5. pc=-4500 

try #15: Target hole radius 1.5 

try #16: Target hole radius .2, Projectile hole radius= .1 
try #17: Etan= 0 in target material model #3 
try #18: Etan= 0, Friction = 0. 

For the large anchor projectile, each file is called a run. Run #1 is the baseline run. Only the 
input data deviations from run #1 are listed. 

run #1: Friction = .l, Target hole radius = .3175, Target material model #3 
run #2: Target material model #5 
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try #7, friction = .1, Constrain target, .7 tgt rad 
dn3d 

Id 1 lp 2 .2 0 .2 35. 

Id 2 lp 2 0 .7 .6 .7 
Id 5 lp 2 2. 35. 2. 5.864 
lap .1 0 -8. 5.864 
Id 14 lp 2 .01 35 2. 35. 

Id 25 lp 2 .02 0 150. 0 
Id 27 lp 2 0 -150. 160 -150. 

Id 22 lp 2 .7 0 .7 -160. 

Id 24 lp 2 15. 0 15. -160. 

Id 28 lp 2 150 0 150 -160. 
si 2 sv fric .1; si 3 tied; 

plane 2 0 0 0 0 1 0 .001 symm 0 0 0 1 0 0 .001 symm 

term .5e-3 plti le-5 

mat 1 3 head 

steel case - 6dac 

ro 7.8e-6 

e 209e6 

pr .3 

sigy 1.2e6 
etan l.e6 
mat 2 3 head 
sandstone 
ro 2.13e-6 
e 13.62e6 
pr .2 

sigy .02e6 
etan .1362e6 
endmat 

part 2 5 14 1 1 2 24 

drag rota 5 0 0 0 0 0 1 90; 

si 2 1 1 2 2 2 2 m velo 0 0 -100e3. 

end 

part 22 27 24 25 2 102 4 
drag rota 5 0 0 0 0 0 1 90; 
si 1 1 1 2 2 1 2 s 
si 1 1 2 2 2 2 3 s 
b 2 1 1222111111 

end 

part 28 25 24 27 2 34 -12 2.0 

drag rota 5 0 0 0 0 0 1 90; 

si 1 1 2 2 2 2 3 m 

b 1 1 2 1 2 1 111111 

end 

end 
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try #8, frict=.l, Larger hole 
dn3d 

Id 1 Ip 2 .2 0 .2 35. 

Id 2 Ip 2 0 .7 .6 .7 
Id 5 Ip 2 2. 35. 2. 5.864 
lap .1 0 -8. 5.864 
Id 14 lp 2 .01 35 2. 35. 

Id 25 lp 2 .02 0 150. 0 
Id 27 lp 2 0 -150. 160 -150. 

Id 22 lp 2 1 0 1 -160. 

Id 24 lp 2 15. 0 15. -160. 

Id 28 lp 2 150 0 150 -160. 
si 2 sv fric . 1; si 3 tied; 

plane 2 0 0 0 0 1 0 .001 symm 0 0 0 1 0 0 .001 symm 

term .5e-3 plti le-5 

mat 1 3 head 

steel case - 6dac 

ro 7.8e-6 

e 209e6 

pr .3 

sigy 1.2e6 
etan l.e6 
mat 2 3 head 
sandstone 
ro 2.13e-6 
e 13.62e6 
pr .2 

sigy .02e6 
etan .1362e6 
endmat 

part 2 5 14 1 1 2 24 
drag rota 5 0 0 0 0 0 1 90; 
si 2 1 1 2 2 2 2 m 
velo 0 0 -100e3. 
end 

part 22 27 24 25 2 102 4 
drag rota 5 0 0 0 0 0 1 90; 
si 1 l 1 2 2 1 2 s 
si 1 1 2 2 2 2 3 s 
b2 1 1 222 111111 
end 

part 28 25 24 27 2 34 -12 2.0 

drag rota 5 0 0 0 0 0 1 90; 

si 1 12 2 2 2 3 m 

bl 12121 111111 

end 

end 
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try #9, frict = . 1, Smaller hole in target 
dn3d 

Id ! Ip 2 .2 0 .2 35. 

Id 2 Ip 2 0 .7 .6 .7 
Id 5 Ip 2 2. 35. 2. 5.864 
lap .1 0 -8. 5.864 
Id 14 lp 2 .01 35 2. 35. 

Id 25 lp 2 .02 0 150. 0 
Id 27 lp 2 0 -150. 160 -150. 

Id 22 lp 2 .51 0 .51 -160. 

Id 24 lp 2 15. 0 15. -160. 

Id 28 lp 2 150 0 150 -160. 
si 2 sv fric .1; si 3 tied; 

plane 2 0 0 0 0 1 0 .001 symm 0 0 0 1 0 0 .001 symm 

term .5e-3 plti le-5 

mat 1 3 head 

steel case - 6dac 

ro 7.8e-6 

e 209e6 

pr .3 

sigy 1.2e6 
etan l.e6 
mat 2 3 head 
sandstone 
ro 2.13e-6 
e 13.62e6 
pr .2 

sigy ,02e6 
etan .1362e6 
endmat 

part 2 5 14 1 1 2 24 
drag rota 5 0 0 0 0 0 1 90; 
si 2 1 1 2 2 2 2 m 
velo 0 0 -100e3. 
end 

part 22 27 24 25 2 102 4 
drag rota 5 0 0 0 0 0 1 90; 
si 1 1 1 2 2 1 2 s 
si 1 1 2 2 2 2 3 s 
b 2 1 1222 111111 
end 

part 28 25 24 27 2 34 -12 2.0 

drag rota 5 0 0 0 0 0 1 90; 

si 1 1 2 2 2 2 3 m 

b 1 12 12 1 111111 

end 

end 


F-5 





try #10, friction = .1, Switch master/slave interface 
dn3d 

Id lip 2 .2 0 .2 35. 

Id 2 lp 2 0 .7 .6 .7 
Id 5 lp 2 2. 35. 2. 5.864 
lap .1 0 -8. 5.864 
Id 14 Ip 2 .01 35 2. 35. 

Id 25 lp 2 .02 0 150. 0 
Id 27 Ip 2 0 -150. 160 -150. 

Id 22 lp2 .7 0 .7 -160. 

Id 24 Ip 2 15. 0 15. -160. 

Id 28 lp 2 150 0 150 -160. 
si 2 sv fric . 1; si 3 tied; 

plane 2 0 0 0 0 1 0 .001 sytnm 0 0 0 1 0 0 .001 symm 

term ,5e-3 plti le-5 

mat 1 3 head 

steel case - 6dac 

ro 7.8e-6 

e 209e6 

pr .3 

sigy 1.2e6 
etan l.e6 
mat 2 3 head 
sandstone 
ro 2.13e-6 
e 13.62e6 
pr .2 

sigy ,02e6 
etan ,1362e6 
endmat 

part 2 5 14 1 1 2 24 

drag rota 5 0 0 0 0 0 1 90; 

si 2 1 1 2 2 2 2 s 

velo 0 0 -100e3. 

end part 22 27 24 25 2 102 4 

drag rota 5 0 0 0 0 0 1 90; 

si 1 1 1 2 2 1 2 m 

si 1 1 2 2 2 2 3 s 

b 2 1 1222111111 

end 

part 28 25 24 27 2 34 -12 2.0 

drag rota 5 0 0 0 0 0 1 90; 

si 1 12 2 2 2 3 m 

b 1 12 12 1 111111 

end 

end 


F-6 


try #11, friction = .05, Constrain target, .7 tgt hole 
dn3d 

Id 1 lp 2 .2 0 .2 35. 

Id 2 lp 2 0 .7 .6 .7 
Id 5 lp 2 2. 35. 2. 5.864 
lap .1 0 -8. 5.864 
Id 14 lp 2 .01 35 2. 35. 

Id 25 lp 2 .02 0 150. 0 
Id 27 lp 2 0 -150. 160 -150. 

Id 22 lp 2 .7 0 .7 -160. 

Id 24 lp 2 15. 0 15. -160. 

Id 28 Ip 2 150 0 150 -160. 
si 2 sv fric .05; si 3 tied; 

plane 2 0 0 0 0 10 .001 symm 0 0 0 1 0 0 .001 symm 

term .5e-3 plti le-5 

mat 1 3 head 

steel case - 6dac 

ro 7.8e-6 

e 209e6 

pr .3 

sigy 1.2e6 
etan l.e6 
mat 2 3 head 
sandstone 
ro 2 13e-6 
e 13.62e6 
pr .2 

sigy .02e6 
etan . 1362e6 
end mat 

part 2 5 14 1 1 2 24 
drag rota 5 0 0 0 0 0 1 90; 
si 2 1 1 2 2 2 2 m 
velo 0 0 -100e3. 
end 

part 22 27 24 25 2 102 4 
drag rota 5 0 0 0 0 0 1 90; 
si 1 1 1 2 2 1 2 s 
si 1 12 2 2 2 3 s 
b 2 1 1 2 2 2 111111 
end 

part 28 25 24 27 2 34 -12 2.0 

drag rota 5 0 0 0 0 0 1 90; 

si 1 12 2 2 2 3 m 

bl 12121 111111 

end 

end 
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try #12, friction = 0, Constrain target, .7 tgt hole 
dn3d 

Id 1 Ip 2 .2 0 .2 35. 

Id 2 lp 2 0 .7 .6 .7 
Id 5 lp 2 2. 35. 2. 5.864 
lap .1 0 -8. 5.864 
Id 14 lp 2 .01 35 2. 35. 

Id 25 lp 2 .02 0 150. 0 
Id 27 lp 2 0 -150. 160 -150. 

Id 22 lp 2 .7 0 .7 -160. 

Id 24 Ip 2 15. 0 15. -160. 

Id 28 lp 2 150 0 150 -160. 
si 2 sv; si 3 tied; 

plane 2 0 0 0 0 1 0 .001 symm 0 0 0 1 0 0 .001 symm 

term .5e-3 plti le-5 

mat 1 3 head 

steel case - 6dac 

ro 7.8e-6 

e 209e6 

pr .3 

sigy 1.2e6 
elan l.e6 
mat 2 3 head 
sandstone 
ro 2.13e-6 
e 13.62e6 
pr .2 

sigy .02e6 
etan . 1362e6 
endmat 

part 2 5 14 1 1 2 24 
drag rota 5 0 0 0 0 0 1 90; 
si 2 1 1 2 2 2 2 m 
velo 0 0 -100e3. 
end 

part 22 27 24 25 2 102 4 
drag rota 5 0 0 0 0 0 1 90; 
si 1 1 1 2 2 1 2 s 
si 1 1 2 2 2 2 3 s 
b 2 1 1 2 2 2 111111 
end 

part 28 25 24 27 2 34 -12 2.0 

drag rota 5 0 0 0 0 0 1 90; 

si 1 1 2 2 2 2 3 m 

b 1 I 2 1 2 1 111111 

end 

end 
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try #13, material model 5 
dn3d 

Id 1 lp 2 .2 0 .2 35. 

Id 2 lp 2 0 .7 .6 .7 
Id 5 lp 2 2. 35. 2. 5.864 
lap .1 0 -8. 5.864 
Id 14 lp 2 .01 35 2. 35. 

Id 25 lp 2 .02 0 150. 0 
Id 27 lp 2 0 -150. 160 -150. 

Id 22 lp 2.7 0.7 -160. 

Id 24 lp 2 15. 0 15. -160. 

Id 28 lp 2 150 0 150 -160. 
si 2 sv fric . 1; si 3 tied; 

plane 2 0 0 0 0 1 0 .001 symm 0 0 0 1 0 0 .001 symm 

term .5e-3 plti !e-5 

mat 1 3 head 

steel case - 6dac 

ro 7.8e-6 

e 209e6 

pr .3 

sigy 1.2e6 
etan l.e6 
endmat 
mat 2 5 head 
sandstone 
ro 2.13e-6 
g 5.68e6 
ku 11.25e6 
aO 1.33 
al 60270 
a2 0 
npts 6 

vs 0 .0025 .01 .04 .05 .055 

p 0 15000 21250 34000 41300 50000 

endmat 

part 2 5 14 1 1 2 24 
drag rota 5 0 0 0 0 0 1 90; 
si 2 1 1 2 2 2 2 m 
velo 0 0 -100e3. 
end 

part 22 27 24 25 2 102 4 
drag rota 5 0 0 0 0 0 1 90; 
si 1 1 1 2 2 1 2 s 
si 1 1 2 2 2 2 3 s 
b 2 1 1 2 22 111111 
end 

part 28 25 24 27 2 34 -12 2.0 

drag rota 5 0 0 0 0 0 1 90; 

si 1 1 2 2 2 2 3 m 

b 1 1 2 1 2 1 111111 

end 

end 
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try #14, material model 5 
dn3d 

Id 1 lp 2 .2 0 .2 35. 

Id 2 lp 2 0 .7 .6 .7 
Id 5 lp 2 2. 35. 2. 5.864 
lap .1 0 -8. 5.864 
Id 14 lp 2 .01 35 2. 35. 

Id 25 lp 2 .02 0 150. 0 
Id 27 lp 2 0 -150. 160 -150. 

Id 22 Ip 2 .7 0 .7 -160. 

Id 24 lp 2 15. 0 15. -160. 

Id 28 lp 2 150 0 150 -160. 
si 2 sv fric .1; si 3 tied; 

plane 2 0 0 0 0 1 0 .001 symm 0 0 0 1 0 0 .001 symm 

term .5e-3 plti le-5 

mat 1 3 head 

steel case - 6dac 

ro 7.8e-6 

e 209e6 

pr .3 

sigy 1.2e6 
etan l.e6 
end mat 
mat 2 5 head 
sandstone 
ro 2.13e-6 
g 5.68e6 
ku 11.25e6 
aO 1.33 
at 60270 
a2 0 

pc -4500 
npts 6 

vs 0 .0025 .01 .04 .05 .055 

p 0 15000 21250 34000 41300 50000 

endmat 

part 2 5 14 1 1 2 24 
drag rota 5 0 0 0 0 0 1 90; 
si 2 1 1 2 2 2 2 m 
velo 0 0 -100e3. 
end 

part 22 27 24 25 2 102 4 
drag rota 5 0 0 0 0 0 1 90; 
si 1 1 1 2 2 1 2 s 
si 1 1 2 2 2 2 3 s 
b 2 1 1 2 2 2 111111 
end 

part 28 25 24 27 2 34 -12 2.0 

drag rota 5 0 0 0 0 0 1 90; 

si 1 12 2 2 2 3 m 

b 1 1 2 1 2 1 111111 

end 

end 
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try #15, frict 3n = .1, Constrain target, 1.5 hole 
dn3d 

Id 1 lp 2 .2 0 .2 35. 

Id 2 lp 2 0 .7 .6 .7 
Id 5 lp 2 2. 35. 2. 5.864 
lap .1 0 -8. 5.864 
Id 14 lp 2 .01 35 2. 35. 

Id 25 lp 2 .02 0 150. 0 
Id 27 lp 2 0 -150. 160 -150. 

Id 22 lp 2 1.5 0 1.5 -160. 

Id 24 lp 2 15. 0 15. -160. 

Id 28 lp 2 150 0 150 -160. 
si 2 sv fric . 1; si 3 tied; 

plane 2 0 0 0 0 1 0 .001 symm 0 0 0 1 0 0 .001 symm 

term .5e-3 plti le-5 

mat 1 3 head 

steel case - 6dac 

ro 7.8e-6 

e 209e6 

pr .3 

sigy 1.2e6 
etan l.e6 
mat 2 3 head 
sandstone 
ro 2.13e-6 
e 13.62e6 
pr .2 

sigy .02e6 
etan . 1362e6 
endmat 

part 2 5 14 1 1 2 24 
drag rota 5 0 0 0 0 0 1 90; 
si 2 1 1 2 2 2 2 m 
velo 0 0 -100e3. 
end 

part 22 27 24 25 2 102 4 
drag rota 5 0 0 0 0 0 1 90; 
si 1 1 1 2 2 1 2 s 
si 1 1 2 2 2 2 3 s 
b 2 1 1 222 111111 
end 

part 28 25 24 27 2 34 -12 2.0 

drag rota 5 0 0 0 0 0 1 90; 

si 1 1 2 2 2 2 3 m 

b 1 12 12 1 111111 

end 

end 
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try #16, friction = .1, Constrain target, .2 tgt hole 
dn3d 

Id 1 Ip 2 .1 0 .1 35. 

Id 2 Ip 2 0 .05 .6 .05 
Id 5 Ip 2 2. 35. 2. 5.864 
lap . 1 0 -8. 5.864 
Id 14 Ip 2 .01 35 2. 35. 

Id 25 lp 2 .02 0 150. 0 
Id 27 lp 2 0 -150. 160 -150. 

Id 22 lp 2 .2 0 .2 -160. 

Id 24 Ip 2 15. 0 15. -160. 

Id 28 lp 2 150 0 150 -160. 
si 2 sv fric . 1; si 3 tied; 

plane 2 0 0 0 0 1 0 .001 symm 0 0 0 1 0 0 .001 symm 

term .5e-3 plti le-5 

mat 1 3 head 

steel case - 6dac 

ro 7.8e-6 

e 209e6 

pr .3 

sigy 1.2e6 
etan t.e6 
mat 2 3 head 
sandstone 
ro 2.13e-6 
e 13.62e6 
pr .2 

sigy ,02e6 
Uan . 1362e6 
endmat 

part 2 5 14 1 1 2 24 
drag ro'a 5 0 0 0 0 0 1 90; 
si 2 1 1 2 2 2 2 m 

velo 0 0 -100e3. 

\ 

end 

part 22 21 24 25 2 102 4 
drag rota 5 0 0 0 0 0 1 90; 
si 1 1 1 2 2 1 2 s 
si 1 1 2 2 2 2 3 s 
b 2 1 1 2 2 2 111111 
end 

part 28 25 24 27 2 34 -12 2.0 

drag rota 5 0 0 0 0 0 1 90; 

si 1 1 2 2 2 2 3 m 

bl 12121 111111 

end 

end 
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try #17, friction = ,1,0 hardning modulus, .7 tgt 
dn3d 

Id 1 Ip 2 .2 0 .2 35. 

Id 2 lp 2 0 .7 .6 .7 
Id 5 lp 2 2. 35. 2. 5.864 
lap .1 0 -8. 5.864 
Id 14 lp 2 .01 35 2. 35. 

Id 23 lp 2 .02 0 150. 0 
Id 27 lp 2 0 -150. 160 -150. 

Id 22 lp 2 .7 0 .7 -160. 

Id 24 lp 2 15. 0 15. -160. 

Id 28 lp 2 150 0 150 -160. 
si 2 sv fric . 1; si 3 tied; 

plane 2 0 0 0 0 1 0 .001 symm 0 0 0 1 0 0 .001 symm 

term .5e-3 plti le-5 

mat 1 3 head 

steel case - 6dac 

ro 7.8e-6 

e 209e6 

pr .3 

sigy 1.2e6 
etan l.e6 
mat 2 3 head 
sandstone 
to 2.13e-6 
e 13.62e6 
pr .2 

sigy .02e6 
etan 0 
end mat 

part 2 5 14 1 1 2 24 
drag rota 5 0 0 0 0 0 1 90; 
si 21 12222m 
velo 0 0 -100e3. 
end 

part 22 27 24 25 2 102 4 
drag rota 5 0 0 0 0 0 1 90; 
si 1 1 1 2 2 1 2 s 
si 1 1 2 2 2 2 3 s 
b 2 1 1 2 2 2 111111 
end 

part 28 25 24 27 2 34 -12 2.0 

drag rota 5 0 0 0 0 0 1 90; 

si 1 1 2 2 2 2 3 m 

b 1 12 12 1 111111 

end 

end 
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try #18, friction = 0 hardening modulus=0 
dn3d 

Id 1 Ip 2 .2 0 .2 35. 

Id 2 lp 2 0 .7 .6 .7 
Id 5 lp 2 2. 35. 2. 5.864 
lap .1 0 -8. 5.864 
Id 14 lp 2 .01 35 2. 35. 

Id 25 lp 2 .02 0 150. 0 
Id 27 lp 2 0 -150. 160 -150. 

Id 22 lp 2 .7 0 .7 -160. 

Id 24 Ip 2 15. 0 15. -160. 

Id 28 lp 2 150 0 150 -160. 
si 2 sv fric 0; si 3 tied; 

plane 2 0 0 0 0 1 0 .001 synun 0 0 0 1 0 0 .001 symm 

term .5e-3 plti le-5 

mat 1 3 head 

steel case - 6dac 

ro 7.8e-6 

e 209e6 

pr .3 

sigy 1.2e6 
etan l.e6 
mat 2 3 head 
sandstone 
ro 2. l3e-6 
e 13.62e6 
pr .2 

sigy ,02e6 
etan 0 
endmat 

part 2 5 14 1 1 2 24 
drag rota 5 0 0 0 0 0 1 90; 
si 2 1 1 2 2 2 2 m 
velo 0 0 -100e3. 
end 

part 22 27 24 25 2 102 4 
drag rota 5 0 0 0 0 0 1 90; 
si 1 1 1 2 2 1 2 s 
si 1 1 2 2 2 2 3 s 
b 2 1 1 222 111111 
end 

part 28 25 24 27 2 34 -12 2.0 

drag rota 5 0 0 0 0 0 1 90; 

si 1 1 2 2 2 2 3 m 

b 1 1 2 1 2 1 111111 

end 

end 
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Run #1, Large Anchor Projectile Baseline 

c friction = . 1 

c target material model = 3 

c target hole radius= .3175 

c units are g, cm, sec. 

dn3d 

c states that a DYNA3D input file will be created 

c next 5 lines describe 1/4 of the projectile profile 

Id 1 lp 2 .254 0 .254 17.78 

Id 2 lp 2 0 .508 1.27 .508 

Id 3 Ip 2 1.27 17.78 1.27 4.211 

lap 0 0 -6.35 4.211 

Id 4 Ip 2 0 17.78 1.27 17.78 

c next 5 lines describe 1/4 of the target profile 

Id 21 lp 2 0. 0. 25.4 0. 

Id 22 lp 2.3175 0.3175 -63.5 
Id 23 lp 2 0. -63.5 25.4 -63.5 
Id 24 lp 2 10.16 0 10.16 -63.5 
Id 26 lp 2 25.4 0. 25.4 -63.5 
c Define the sliding interfaces 
si 1 sv fric .1; si 2 tied; 

c Sliding interface 1 is sliding with voids and has . 1 coef. of friction 

c Sliding interface 2 is a tied boundary 

plane 2 0 0 0 0 1 0 .001 symm 0 0 0 1 0 0 .001 symm 

c describes the planes of symmetry 

c The termination time and plot interval are set. 

term 5e-3 plti le-5 

c Start define material #1 

mat 1 3 head 

4340 steel projectile 

c assign material model #3 and puts in a heading statement 
c density = ro = g/cm A 3 
ro 7.8 

c Young’s modulus = e = g/(cm*s*2) 
e 2.00el2 

c Poison’s ratio = pr 
pr .3 

c Yield stress = sigy = g/(cm*s'2) 
sigy 17.9e9 

c Tangential modulus = etan = g/(cm*s*2) 

etan 10.0e9 

mat 2 3 head 

sandstone 

ro 2.16 

e 44e9 

pr .3 

sigy .234e9 
etan .440e9 
end mat 

part 2 3 4 1 1 2 14 
drag rota 5 0 0 0 0 0 1 90; 
si 2 1 1 2 2 2 1 m 

c initial velocity of the projectile in cm/sec. 
velo 0 0 -10000. 
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end 

part 22 23 24 21 2 100 4 
drag rota 5 0 0 0 0 0 1 90; 
si 1 1 1 2 2 1 1 s 
si 1 1 2 2 2 2 2 s 
b 2 1 1 2 2 2 111111 
end 

part 26 21 24 23 2 25 -3 2.0 

c defines a four sided region, lines 26, 21, 24, 23, assigns material 2 

c 34 elements between lines 21 & 23, 12 elements between lines 24 & 26 

c -12 (- sign) & the 2 defines the ratio of the first to the last 

c segment between lines 24 and 26. 

drag rota 5 0 0 0 0 0 1 90; 

si 1 1 2 2 2 2 2 m 

b 1 12 12 1 111111 

end 

end 
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Run Wl, Large Anchor Projectile Baseline 
c target material model = S 
c units are g, cm, sec. 
dn3d 

c states that a DYNA3D input file will be created 

c next 5 lines describe 1/4 of the projectile profile 

Id 1 Ip 2 .254 0 .254 17.78 

Id 2 lp 2 0 .508 1.27 .508 

Id 3 Ip 2 1.27 17.78 1.27 4.211 

lap 0 0 -6.35 4.211 

Id 4 lp 2 0 17.78 1.27 17.78 

c next 5 lines describe 1/4 of the target profile 

Id 21 lp 2 0. 0. 25.4 0. 

Id 22 lp 2 .3175 0 .3175 -63.5 
Id 23 lp 2 0. -63.5 25.4 -63.5 
Id 24 lp 2 10.16 0 10.16 -63.5 
Id 26 lp 2 25.4 0. 25.4 -63.5 

c Define the sliding interfaces si 1 sv fric .1; si 2 tied; 

c Sliding interface 1 is sliding with voids and has . 1 coef. of friction 

c Sliding interface 2 is a tied boundary 

plane 2 0 0 0 0 1 0 .001 symm 0 0 0 1 0 0 .001 symm 

c describes the planes of symmetry 

c The termination time and plot interval are set. 

term 5e-3 plti le-5 

c Start define material #1 

mat 1 3 head 

4340 steel projectile 

c assign material model #3 and puts in a heading statement 
c density = ro = g/cm A 3 
ro 7.8 

c Young’s modulus = e = g/(cm*s A 2) 
e 2.00el2 

c Poison’s ratio = pr 
pr .3 

c Yield stress = sigy = g/(cm*s A 2) 
sigy 17.9e9 

c Tangential modulus = etan = g/(cm*s A 2) 
etan 10.0e9 
mat 2 5 head 

sandstone, material model ft 5 
ro 2.16 

c Shear modulus = g = g/(cm*s A 2) 
g 8.44el0 

c Bulk unloading modulus = ku = g/(cm' ,, s A 2) 
ku 11.25el0 

c a0-2 are constants for the Yield function 
aO ,00133e7 
al 60.27e7 
a2 0 

c pc = pressure cutoff 
pc 0 
npts 6 

c npts specifies 6 points on the p-vs curve 
c vs = volumetric strain 


F-17 





vs 0 .0025 .01 .04 .05 .055 
c p = pressure 

p 0 15e7 21.25e7 34e7 41.3e7 50e7 
endmat 

part 2 3 4 1 1 2 14 
drag rota 5 0 0 0 0 0 1 90; 
si 2 1 1 2 2 2 1 m 

c initial velocity of the projectile in cm/sec. 

velo 0 0 -10000. 

end 

part 22 23 24 21 2 100 4 
drag rota 5 0 0 0 0 0 1 90; 
si 1 1 1 2 2 1 1 s 
si 1 1 2 2 2 2 2 s 
b 2 1 1 2 2 2 111111 
end 

part 26 21 24 23 2 25 -3 2.0 

c defines a four sided region, lines 26, 21, 24, 23, assigns material 2 

c 34 elements between lines 21 & 23, 12 elements between lines 24 & 26 

c -12 (- sign) & the 2 defines the ratio of the First to the last 

c segment between lines 24 and 26. 

drag rota 5 0 0 0 0 0 1 90; 

si 1 1 2 2 2 2 2 m 

b 1 1 2 1 2 1 111111 

end 

end 
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